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INTRODUCTION 


At  the  forefront  of  modern  research  in  structural  mechanics  is  the 
problem  of  dealing  with  time-dependent,  nonlinear  material  models.  The 
stimulus  for  this  interest  is  due  to  the  increased  use  of  new  materials 
in  structural  systems,  such  as  plastics,  glues,  resins,  biological 
tissues,  fuel  propellants,  soils,  sea-ice,  and  high  temperature  metals. 

Although  constitutive  theories  for  time-dependent,  nonlinear  models 
are  still  in  their  infancy,  the  ability  to  solve  boundary  value  problems 
embracing  these  theories  is  well  in  hand.  In  particular,  the  well-known 
finite  element  method  is  capable  of  incorporating  the  most  elegant 
constitutive  models  with  relative  ease  and  sufficient  accuracy  for 
engineering  applications.  This  is  a complete  turn  of  events  from 
classical  mechanics,  wherein  constitutive  theories  were  necessarily 
simplified  to  afford  tractable  solutions  of  boundary. value  problems. 

Clearly,  the  time  has  come  to  re-focus  attention  on  material 
modeling  to  develop  constitutive  theories  that  are  representative  of 
real  material  behavior  and  are  on  a par  with  modern  analytical  methods. 


BACKGROUND 

Time-dependent  and  nonlinear  constitutive  models  can  generally  be 
identified  as  belonging  to  one  of  three  classifications:  creep,  non- 
linear viscoelastic,  or  viscoplastic.  The  first  two  classifications  are 
fairly  well  defined  in  the  literature  [1,2]  and  will  not  be  pursued 
here.  Suffice  to  say,  these  models  have  met  with  limited  success  in 
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some  engineering  applications,  but  they  are  not  able  to  replicate  many 
of  the  observed  phenomena  of  real  materials,  such  as  the  Bauschinger 
effect,  anisotropic  deformation,  and  cyclic  phenomena  [3].  f 

y j 

As  the  name  suggests,  a viscoplastic  model  combines  the  features  of  \ 
viscous  time-dependency  with  plasticity  theory,  and  in  its  general  form  } 
is  capable  of  faithfully  representing  observed  material  behavior,  such  | 

as  a rate-dependent  yield  strength,  Bauschinger  effect,  and  loading  and  | 

* 

i 

unloading  phenomena.  These  concepts  will  be  amplified  throughout  the  ^ 

f 

text.  4 

s 

The  genera]  theory  of  viscoplasticity  is  well  summarized  by  Perzyna  js 
[A],  wherein  he  concludes  further  theoretical  development  is  necessary.  | 

i 

However,  two  important  subsets  of  the  general  theory  appear  theoreti-  § 

ealLy  sound  and  tractable.  The  first  will  be  called  "clastic-viscoplastic" 

i 

and  the  second,  "viscoelastic-plastic."  The  former  is  characterized  by 

an  elastic  regLon  within  the  yield  surface  and  a time-dependent  yield  f 

I 

function  whose  domain  is  not  restricted  to  the  yield  surface.  ZlenkiuwLczi 

I 

[5]  has  developed  and  demonstrated  this  model  in  several  finite  element  J 

f 

applications.  ‘ 

The  "viscoelastic-plastic"  model  is  characterized  by  linear  visco-  ? 

elastic  theory  within  the  yield  surface  and  the  combined  influence  of  | 

viscoelasticity  and  plasticity  on  the  vipld  surface.  The  yield  function  | 

f 

is  restricted  to  the  domain  within  and  on  the  yield  surface  as  in  classi-  « 
cal  plasticity.  This  model  is  developed  in  detail  in  this  report.  | 


OBJECTIVE  AND  SCOPE 

The  objective  of  this  study  is  to  present  an  incremental  constitu- 
tive development  for  the  viscoelastic-plastic  model,  and  furthermore,  to 
incorporate  the  model  into  a finite  element  formulation.  To  achieve 
this  objective,  the  classical  theories  of  plasticity  and  linear  visco- 
elasticity are  individually  developed  in  preparation  for  the  combined 
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viscoelastic-plastic  development.  In  themselves,  the  individual  treat- 
ments contain  many  generalizations  and  unifications  of  plasticity  and 
viscoelasticity  heretofore  scattered  throughout  the  literature.  Conse- 
quently, this  work  serves  as  a reference  to  plasticity  and  visco- 
elasticity in  addition  to  viscoplasticity. 

It  is  well  recognized  that  any  constitutive  theory  is  useful  only 
if  it  can  be  successfully  incorporated  into  the  field  equations  of 
boundary  value  problems.  Accordingly,  each  of  the  models  presented  in 
this  work  is  accompanied  with  a general  finite  element  formulation  and 
suggested  solution  algorithms.  Since  the  finite  element  formulation  is 
common  to  all  material  models,  it  will  be  presented  at  the  outset  of 
this  study,  thereby  providing  some  insight  into  the  desired  form  of  the 
constitutive  relationships. 
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Chapter  1 


# 


NONLINEAR  FINITE  ELEMENT  FORMULATION  OF  EQUILIBRIUM  EQUATIONS 


The  purpose  of  this  section  is  to  formulate  a finite  element 
methodology  capable  of  handling  all  of  the  nonlinear  constitutive  models 
presented  later  in  the  text.  At  first  thought,  this  may  appear  as  put- 
ting the  "cart  before  the  horse."  However,  this  ordering  of  the  presen- 
tation provides  the  reader  with  an  informative  preview  of  the  assumptions 
common  to  all  constitutive  theories,  as  well  as  demonstrating  how  the 
various  constitutive  laws  are  utilized  to  solve  boundary  value  problems. 

Perhaps  the  most  powenui  piiysical  law  in  analytical  mechanics  is 
the  principle  of  virtual  work.  By  evoking  this  law  under  isothermal 
conditions,  the  resulting  statement  is  completely  general  and  valid  for 
both  geometric  and  material  nonlinearities.  To  wit,  virtual  work  may  be 
stated  as:  Given  a deformed  body  in  equilibrium  under  a set  of  external 
loads,  and  subjected  to  any  small  virtual  displacement  compatible  with 
the  constraints  of  the  deformed  body;  then,  the  virtual  work  of  the 

internal  forces  is  equal  to  the  virtual  work  of  the  external  loads. 

Furthermore,  by  d'Alembert's  principle,  inertia  loads  may  be  treated  as 

equivalent  external  loads,  thereby  extending  the  virtual  work  concept  to 
dynamic  problems.  Tn  equation  form  this  statement  is: 


J'  o dv  = J"  6u^  t ds  + y 6u^(f  - ii 


u ) p dv 


where 


a = stress  vector  with  respect  to;?flef ormed  bodv  (i.e., 
Cauchy  stress) 


O-l)  i 
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= infinitesimal  strain  vector  with  respect  to  deformed 
body 

= displacement  vector  measured  from  some  fixed  inertial 
Cartesian  reference  system  (u  denotes  acceleration) 

= surface  traction  vector  on  deformed  body 

= current  mass  density  per  unit  volume 

= body  force  vector  per  unit  mass 

= current  volume  of  body 

= current  surface  of  body 

=)  variational  quantity 

=)  vector  quantity 

=)  transpose  of  a vector 

With  regard  to  large  deformations,  Equation  1-1  poses  a rather 
difficult  nonlinear  problem  since  all  quantities  are  referenced  to  the 
unknown  deformed  configuration.  The  above  formulation  of  the  problem  is 
often  called  the  "Eulerian  description,"  or  "moving  coordinate  approach." 
Alternatively,  it  is  possible  to  recast  the  formulation  so  that  all  quan- 
tities are  measured  with  respect  to  the  original  or  initial  configuration 
this  is  known  as  the  "Lagrangian  description."  By  this  description,  all 
derivatives  and  integrals  are  easily  computed  with  respect  to  the  initial 
fixed  coordinate  system.  However,  the  stress  and  strain  vectors  take  on 
new  definitions  consistent  with  their  reference  to  the  original  configu- 
ration. Namely,  they  become  the  second  Piola-Kirchof f stress  vector  and 
the  Creen-l.agrange  strain  vector  [6],  both  of  which  are  nonlinear  with 
respect  to  the  deformation  gradient. 

Eulerian  and  Lagrangian  descriptions  represent  the  end  points  on  a 
spectrum  of  possible  formulations  for  large  deformation  analysis.  Many 
innovative  combinations  and  approximations  of  these  formulations  have 


6(  ) 
( ) 
( )T 
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been  presented  in  the  literature  [7,8,9],  However,  the  most  significant 
point  to  be  kept  in  mind,  with  regard  to  both  geometric  and  material 
nonlinearity,  is  that  the  constitutive  law  must  be  defined  in  a manner 
consistent  with  the  formulation.  For  example,  a Lagrangian  formulation 
requires  that  the  law  relate  second  Piola-Kirchof f stress  to  Green- 
Lagrange  strain. 

For  the  remainder  of  this  development,  infinitesimal  strain  theory 
will  be  assumed  and  inertia  loads  will  be  neglected  so  that  attention 
can  be  focused  on  material  nonlinearity.  Nonetheless,  the  constitutive 
models  presented  herein  are  equally  valid  for  large  deformations  and 
dynamic  loadings,  provided  they  are  used  in  a consistent  fashion  as 
noted  above. 

In  accordance  with  infinitesimal  strain  theory,  the  deformed  con- 
figuration is  assumed  to  differ  infinitesimally  from  the  undeformed  * • 
configuration;  consequently,  Equation  1-1  Is  assumed  valid  for  all  1 
quantities  measured  with  respect  to  the  undeformed  configuration. 
Accordingly,  the  stress  and  strain  vectors  have  the  classical  defi- 
nition. 

In  addition  to  Equation  1-1,  the  complete  formulation  of  a boundary 
value  problem  requires  the  specification  of  (1)  the  strain-displacement 
relationship,  (2)  constitutive  law  (i.e.,  the  stress-strain  relationship), 
and  (3)  the  boundary  conditions.  For  the  present,  no  restrictions  on 
the  form  of  the  constitutive  law  will  be  assumed;  however,  the  strain- 
displacement  relationship  will  be  assumed  linear  in  accordance  with 
infinitesimal  strain  theory,  and  the  boundary  conditions  will  be  intro- 
duced via  the  finite  element  method. 

Thus,  the  strain-displacement  relationship  takes  the  form: 


e = Ou 


(1-2) 
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where  0 is  a linear  operator  on  u.  For  example,  in  a Cartesian  coordi- 

M ~ 

nate  system,  0 has  the  form  such  that  the  normal  strain  components  are 

m 

given  by  e, . = 3u^/3X^  (no  sum  on  i) , and  the  shear  strain  components 

are  given  by  2c..  = y.  . = 3u./3X.  + 3u./3X.. 

ij  ij  i J J i 

At  this  juncture,  the  fundamental  finite  element  approximation  is 
introduced.  The  body  is  subdivided  into  a discrete  set  of  elements 
whereby  the  integrations  in  Equation  1-1  apply  to  each  element.’  Within 
each  element,  the  displacement  vector  is  approximated  by  an  assumed  • 
spatial  interpolation  function.  Associated  with  the  interpolation 
function  are  a set  of  unknown  nodal  displacements  located  on  the  surface 
of  each  element.  Adjacent  elements  share  common  nodes,  thereby  pro- 
viding nodal,  point  compatibility  between  the  elements.  If  the  interpo- 
lation function  is  admissible  [101,  it  can  be  shown  that  the  finite 
clement  approximation  will  provide  the  "exact"  solution  to  the  boundary.-' 
value  problem  with  a sufficient  number  of  elements.  Symbolically,  tlTe 
finite  element  approximation  Is  written  ns: 


= h 


A 

U 


0-3) 


where  h = matrix  of  admissible  interpolation  functions  dependent 
only  on  space 

u = vector  of  nodal  displacements  dependent  only  on  time 
(or  load  step) 

Inserting  Equation  1-3  into  Equation  1-2,  the  strain-displacement 
relationship  within  each  element  is  given  as: 


- « A 
15  u 


(1-4) 
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where  B equals  Oh,  which  is  the  stra in-to-nodal  displacement  matrix 

•a  ^ 

dependent  only  on  spatial  coordinates.  Returning  to  the  statement  of 
virtual  work,  and  introducing  the  finite  element  approximations  imbedded 
in  Equations  1-3  and  1-4,  Equation  1-1  can  be  written  as  (neglecting 
inertia) : 


f p dv  + 


0-5) 


Here  the  symbol  'i  denotes  summation  over  the  elements  with  the  special 
understanding  that  the  contributions  of  each  element  are  properly  assigned 
to  the  corresponding  location  of  the  nodal  displacement  vector. 

Since  each  variation  of  u is  independent,  Equation  1-5  represents  a 
set  of  simultaneous  equations.  In  preparation  for  the  incremental  con- 
stitutive relationships  to  be  introduced  next,  and  without  loss  of 
generality.  Equation  1-3  can  be  written  in  incremental  form  by  replacing 
o,  t,  and  f by  3a,  At  and  At,  respectively,  to  give  the  set  of  simul- 
taneous equations: 


s / 


BT  Ao  ch 


AP 


(1-6) 


where 


(1-7) 


In  a similar  manner,  the  strain  displacement  relationship  can  be 
written  in  incremental  form  as: 


8 


Ae 


(1-8) 


= B Au 

•m  ~ 


Thus  far  in  the  discussion  no  assumptions  have  been  made  with 
regard  to  the  form  of  the  constitutive  relationship.  Therefore,  Equa- 
tions 1-6,  1-7,  and  1-8,  which  are  completely  general  with  regard  to 
material  laws,  provide  the  proper  starting  point  for  incorporating 
material  models  into  a finite  element  formulation.  For  present  pur- 
poses, no  special  restrictions  will  be  imposed  on  the  form  of  the  con- 
stitutive relationships.  However,  for  the  sake  of  clarity  and  ease  of 
presentation,  the  constitutive  relationship  will  be  symbolically  denoted 
as : 


Ac  = D At 
/w  * n <— 


(1-9) 


e,  represents  a stress-9train  relationship  that  may  depend  on 

total  stress,  total  strain,  and  history  of  loading.  Accordingly,  D mav 

«»» n 

be  construed  as  an  operator  matrix  rather  than  a simple  matrix. 

The  motivation  for  introducing  the  incremental  quantities  Ao  and 
Ar.  into  Equation  1-9  is  to  facilitate  linearization  of  D . For  example, 

* II 

if  sufficiently  small  time  steps  (or  load  steps)  are  prescribed,  it  may 

be  sufficientlv  accurate  to  determine  D based  on  the  stress-strain 

~n 

state  at  the  beginning  of  the  load  step  (tangent  modulus  approach).  On 
the  other  hand,  if  larger  time  steps  are  prescribed  or  the  material 

character  is  highly  nonlinear,  it  mav  he  neeessarv  to  determine  1)  based 

»n 

on  the  average  stress-strain  state  over  the  interval,  thereby  r''  airing 
iteration  within  the  time  step  (modified  tangent  or  chord  modulus 
approach) . 

In  either  approach,  the  nature  and  the  treatment  of  D are  dis- 
cussed individually  for  each  constitutive  model  presented  in  later 

sections  of  this  report.  For  now,  it  is  simple  assumed  that  D exists 

*n 

and  may  or  may  not  be  dependent  upon  the  current  time  step. 


9 


With  the  above  understanding,  two  basic  solution  strategies  for 
incorporating  the  incremental  constitutive  law  into  the  finite  element 
formulation  are  (1)  the  "tangent  stiffness  method,"  and  (2)  the  "initial 
strain  method."  The  names  of  both  methods  are  misleading;  however,  the 
use  of  this  terminology  Is  widespread  in  engineering  literature  [11,12] 
and  will  be  retained  here.  Broadly  speaking,  the  tangent  stiffness 
method  implies  the  constitutive  relationship  is  incorporated  directly 
into  the  global  stiffness  matrix,  whereas  the  initial  strain  method 
incorporates  the  linear  portion  of  the  constitutive  relationship  into 
the  global  stiffness  matrix,  treating  nonlinear  terms  as  load  vectors. 
Both  methods  are  described  below. 


TANGENT  STIFFNESS  METHOD 


The  tangent  stiffness  method  is  a straightforward,  brute  force 
technique  that  requires  no  further  assumptions  than  already  presented. 

Formally,  Equations  1-8  and  1-9  are  combined  to  give  Ao  = D B Au,  and 

% * n • ^ 

this  result  is  inserted  into  the  global  equilibrium  equations,  Equation 
1-6,  to  give: 


IS1  Ai! 


(1-10) 


where 


E / 


B D B dv 

» % n w 


(1-id 


In  the  above  K,  is  the  global  stiffness  matrix  relating  displace- 

mcnt  itn  etnents  to  load  increments.  Since  D is  dependent  on  the 

•.n 

stress-strain  state,  KT  changes  accordingly. 

Tlie  solution  algorithm  embodies  the  ml  lowing  steps,  wherein  it  is 
assumed  the  system  is  in  equilibrium  at  time  step  "i”  with  known  response.' 
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Oj,  c ^ , and  G The  objective  is  to  find  the  response  increments  Au, 
Ae,  and  Ao  over  the  interval  i to  i + 1. 

1.  Calculate  load  increment  AP  (Equation  1-7). 

2.  Assemble  global  stiffness  K_,  based  on  value  of  D at  beginning 

* i **n 

of  step  (Equation  1-11) . 

3.  Solve  K Au  = AP , for  Au. 

w L **  ~ ~ 

4.  Determine  At  and  Ac  from  Equations  1-8  and  1-9. 

5.  If  desired,  re-evaluate  based  on  average  values  of  stress- 
strain  over  the  interval,  and  return  to  Step  2 to  iterate 
within  time  step. 

6.  Add  incremental  responses  to  total  responses,  advance  the  time 
step,  and  return  to  Step  1. 

In  the  above  algorithm  the  global  matrix  K.r  must  be  assembled  and 
triangular ized  for  each  load  step  and  any  iteration  within  the  load 
step.  For  large  problems  this  method  may  be  prohibitively  costly. 


INITIAL  STRAIN  METHOD 

The  initial  strain  method  generally  provides  a much  more  efficient 
algorithm;  however,  two  additional  assumptions  are  required  with  regard 
to  the  material  behavior.  First,  it  is  assumed  the  strain  increment  can 
be  decomposed  into  Linear  and  nonlinear  contributions.  That  is,  the 
total  strain  increment  is  given  by: 


(1-12) 
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where  Ae^  is  the  linear  strain,  and  Ac^  represents  a combined  lumping  of 
all  nonlinear  strains.  Secondly,  it  is  assumed  the  total  stress  incre- 
ment is  persistently  related  to  the  linear  strain  increment  by: 


lo  = D Ac 

— « 1.  ~L 


(1-13) 


where  D is  a linear  material  matrix.  Combining  Equations  1-12  and  1-13 

M L 

gives: 


Ao  = D (Ae  - At) 


(l-U) 


Recalling  the  general  constitutive  assumption,  Ao  = D At,  and  using 
Equation  1-14,  a relationship  between  total  and  nonlinear  strain  incre- 
ments is  obtained  as: 


1),  At 

*.L  ~n 


(Dn  " D1  )A£ 

w n % L«  /v 


(1-13) 


To  make  use  of  the  above  in  the  initial  strain  formulation,  Ao  of 
Equation  1-6  is  replaced  by  Ao  = (At  - Ac  ) , and  the  nonlinear  terms 
are  brought  to  the  right-hand  side  to  give: 


M>  + AF 


(1-16) 


where 


S/."T  h. 


(1-17) 
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AF 

-n 


(1-18) 


- tzJ 

v 


B D,  At  dv 
m »L  ~n 


Equation  1-16  is  the  global  equilibrium  equation  for  the  initial 
strain  formulation.  is  the  constant  elastic  stiffness  matrix,  and 

AF^  represents  a nonlinear  force  vector  composed  of  nonelastic  strain 
increments . 

Since  AF  is  an  unknown  vector,  it  is  treated  iterativelv.  To  this 
~n 

end,  it  is  convenient  to  use  Equation  1-15  and  the  strain  displacement 
relationship  to  replace  D,  Ac  by  (D,  - D )B  u.  That  is,  AF  mav  be 
written  as: 

AF  = / BT(D,  - D )B  dv  Au  (1-19) 

~n  I jLJ  J - -L  «n  » I ~ 


Thus,  the  basic  idea  is  to  solve  Equation  1-16  for  Au  and  then  redefine 
AF^  accordingly  until  convergence  is  achieved. 

More  formally,  the  basic  steps  of  the  initial  strain  algorithm  are 
given  below.  It  is  assumed  the  system  is  in  equilibrium  at  time  step  i. 
The  objective  is  to  determine  the  response  increments  Au^,  etc.,  for  the 
next  time  step  i + 1. 


1.  Calculate  load  increment  .AT  . . 

- i 

2.  Estimate  AF  = AF 

~n  , ~n  . , 


3. 

So  lve 

(back 

substitute) : VL 

4 A 

:)U  . 

= AP  . + AF 

* L 

-1 

~ l ~n  ^ 

- D )B  dv"!  Au. 

4 , 

Recal 

c u 1 a t e 

. ..  - -k  r f 

B1  0), 

5.  Compare  successive  estimates  of  AFn  or  Au ^ . If  converged,  go 
to  Step  6.  Otherwise,  repeat  iteration  loop  3-4-5. 

6.  Advance  load  step,  and  return  to  Step  1. 

The  most  significant  feature  of  the  initial  strain  algorithm  is 

that  need  only  be  triangular ized  once  at  the  outset  of  the  calculations, 

thereby  allowing  rapid  solutions  by  back,  substitution.  Furthermore, 

within  the  iteration  loop,  the  matrix  D need  not  be  assumed  constant 

n 

but  can  be  simultaneously  modified  in  accordance  with  Au. 

The  initial  strain  algorithm  will  vary  slightly  for  different 
constitutive  theories  (e.g.,  plasticity,  viscoelasticity,  and  visco- 
plasticity). Nonetheless,  the  basic  features  remain  the  same. 


SUMMARY 


The  intent  of  this  section  was  to  provide  the  reader  with  the 

framework  for  implementing  the  constitutive  theories  to  be  discussed. 

Moreover,  it  is  hoped  that  an  appreciation  for  some  common  assumptions 

inherent  in  nonlinear  constitutive  models  was  achieved.  To  wit,  an 

incremental  constitutive  equation  of  the  form  Ao  = D Ac  is  the  desired 

~ »n  ~ 

objective. 


14 


ii,L'  aii  ' -"•iit.b.iiL ..JtaJic  i . wj,..,. .uih«ul:;i  wL..,* iIilmhimmIiLmL i ^ ^ 


Chapter  2 


PLASTICITY  THEORY 


INTRODUCTION 


The  original  intention  of  this  section  was  to  provide  a review  of 
classical  plasticity  theory  for  the  purpose  of  establishing  notation  and 
concepts  to  be  used  later  in  the  presentation  of  viscoplasticity, 
However,  upon  reviewing  plasticity  literature,  no  single  presentation 
was  found  that  could  be  termed  'classical'  in  the  sense  that  it  was  both 
complete  and  conceptually  instructive.  For  example,  introductory  texts 
[11,14]  deal  almost  exclusively  with  a particular  form  of  plasticity, 
e.g.,  Prandtl-Reuss  equations  with  isotropic  hardening.  Thus,  their 
treatment  is  incomplete.  On  the  other  end  of  the  scale,  the  presenta- 
tion offered  by  Nayak  and  Zienkiewicz  [15]  is  a unified,  complete 
plasticity  theory  formulation;  however,  it  offers  lit.  , in  the  way  of 
conceptual  insights.  Furthermore,  nowhere  in  the  literature  is  there  a 
general  treatment  of  the  "universal  hardening  law"  which  combines  kine- 
matic and  isotropic  hardening  into  r.  unified  theory.  This  concept  was 
originally  introduced  by  Hodge  [luj,  and  further  discussed  by  Goel  and 
Malvern  117],  but  it  has  not  yet  been  treated  in  a unified  manner. 

This  presentation  is  an  attempt  at  a unified  plasticity  theory  with 
universal  hardening.  Moreover,  an  original  one-dimensional  model  is 
presented  that  offers  an  insight  into  the  nature  of  isotropic,  kine- 
matic, and  universal  hardening.  Accordingly,  it  is  hoped  the  develop- 
ments herein  will  offer  more  than  a simple  review. 
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The  final  objective  of  this  plasticity  presentation  is  to  obtain  a 
constitutive  relationship  compatible  with  the  finite  element  formulations 
previously  developed.  In  particular,  the  following  relationship  is 
desired : 


Ao  = D Ac 

„ep  ~ 


(2-1) 


where 


Ao  = stress  increment  vector 


Ac  = strain  increment  vector 


D = incremental  elasti. -plastic  constitutive  matrix 
«ep 


To  achieve  this  objective,  a rather  detailed  development  of  plas- 
ticity theory  is  offered  which  ultimately  produces  the  desired  matrix 

D . Lastlv,  the  initial  strain  finite  element  formulation  for  plasti- 
»ep 

city  is  presented. 


PLASTICITY  ASSUMPTIONS  AND  CONDITIONS 


Common  to  all  plasticity  developments  is  the  assumption  that  total 
strain  increment  can  be  decomposed  into  elastic  and  plastic  components, 
and  further,  it  is  assumed  the  stress  increment  is  persistently  related 
to  the  elastic  strain  increment  through  generalized  Hook's  law.  That 


where  Ac 
~e 


Aa  = D Ae 
*e  ~e 


= elastic  strain  increment 


= elastic  constitutive  matrix  (Hook’s  law) 


(2-3) 


The  above  equations  are  universally  assumed  in  all  plasticity 
formulations  known  to  the  author,  and  indeed,  are  identical  to  the 
assumptions  employed  in  the  finite  element  initial  strain  formulations. 
In  the  course  of  this  work  some  insights  into  the  motivation  of  these 
assumptions  will  be  given. 

It  is  important  to  note  that  the  developments  presented  herein 
assume  plastic  hardening.  The  special  case  of  no  hardening  (perfect- 
plasticity)  is  easily  dealt  with  once  the  general  equations  have  been 
developed  on  the  presumption  of  hardening. 

In  addition  to  the  preceding  assumptions,  plasticity  theory  is 
built  upon  four  basic  conditions.  These  conditions  are  formally  defined 
below  and  will  be  clarified  in  the  ensuing  discussion. 

1*  Yield  Condition.  A scalar  function  signifying  when  plastic 
yielding  will  occur.  It  is  composed  of  a positive  valued  "loading 
function”  and  a "yield  parameter"  that  is  never  less  than  the  value  of 
the  loading  function.  Yielding  can  only  occur  when  the  loading  function 
is  equal  to  the  yield  parameter.  The  states  of  stress  satisfying  this 
condition  form  a hypersurface  in  stress  space,  terme.  .he  yield  surface. 

2.  Surface  Hardening.  A law  that  "tracks"  or  measures  the  move- 
ment of  the  yield  surface  during  plastic  yielding.  Traditional!  , yield 
surface  movement  is  restricted  to  uniform  expansion  (isotropic  hardening) 
or  "rigid  body"  translation  (kinematic  hardening).  The  former  is 
denoted  by  an  increase  in  the  yield  parameter,  while  the  latter  is 
denoted  by  a translation  in  stress  space  of  the  loading  function.  The 
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linear  combination  of  kinematic  and  isotropic  hardening  constitutes 
"universal  hardening."  (Note,  perfect-plasticity  implies  no  hardening, 
l.e.,  yield  surface  never  changes). 

3.  Flow  Rule.  An  expression  relating  increments  of  plastic  strain 
to  increments  of  stress.  The  magnitude  of  the  incremental  plastic 
strain  vector  is  dependent  upon  the  magnitude  of  the  stress  increment 
normal  to  the  yield  surface  and  a modulus  known  as  the  "hardening 
coefficient."  The  direction  of  the  incremental  plastic  strain  vector 
may  he  assumed  normal  to  the  yield  surface  (associative  law),  or  normal 
to  some  other  hypersurface  (nonassociative) . 


U.  Hardening  Rule.  An  expression  or  set  of  data  points  that  pro- 
vides the  value  of  the  "hardening  coefficient"  as  a function  of  some 
plastic  deformation  measure,  such  as  plastic  work  or  effective  plastic 
strain.  (Note,  this  is  distinct  from  the  concept  of  "surface  hardening" 
described  ahove  because  one  can  specify  a hardening  rule  Independent  of 
'.Tie  form  of  surface  hardening.) 


In  order  to  amplify  and  clarify  these  four  concepts,  a series  of 
one-dimensional  conceptual  models  are  introduced  demonstrating  the 
nature  of  kinematic,  isotropic,  and  universal  surface  hardening.  • For 
each  of  the  models,  the  associated  one-dimensional  yield  condition  is 
generalized  for  multidimensional  stress  states.  Also,  the  corresponding 
rules  for  "tracking"  the  current  location  and  size  of  the  yield  surface 
are  presented.  Although  the  universal  model  embraces  both  the  kinematic 
and  isotropic  models,  it  is  instructive  to  study  the  kinematic  and 
isotropic  models  individually  to  better  appreciate  the  universal  model. 

The  concept  of  the  "flow  rule"  and  the  associated  "hardening  rule" 
is  presented  after  the  generalized  yield  condition  with  universal  hard- 
ening. Lastly,  having  established  the  four  basic  concepts  of  plasticity, 
expressions  for  D are  developed  for  the  solution  algorithms.  Supple- 
mentary discussions  on  loading  functions  and  nonassociative  flow  ruleg 
are  given  in  Appendix  A.  ■ ' 


PLASTICITY  CONCEPTS  WITH  ONE-DIMENSIONAL  MODELS 


Plastic  deformation  is  often  called  frictional  deformation  because 
of  the  similarity  between  plasticity  and  the  classical  Coulomb  frictional 
hypothesis  for  sliding  bodies.  To  wit,  plastic  deformation  (sliding 
movement)  only  occurs  when  the  yield  stress  (frictional  resistance)  is  * 
exceeded  by  active  loads.  Moreover,  during  plastic  deformation,  the 
plastic  work  (dissipated  energy)  is  independent  of  the  rate  of  deforma- 
tion (rate  of  sliding)  and  only  dependent  on  the  deformation  path  (sli- 
d ing  path) . 

This  analogy  between  frictional  theory  and  plasticity  theory  can  be 
exploited  to  produce  conceptual  models  that  provide  a keen  insight  into 
the  behavior  of  elastic-plastic  materials. 


Kinematic  Model 


Figure  la  represents  a one-dimensional  elastic-plastic  model  with 
kinematic  hardening.  The  sliding  block  obeys  the  Coulomb  friction 
hypothesis  in  that  the  frictional  resistance,  o^.,  is  a passive  resis- 
tance, i.e.,  equal  but  opposite  to  the  unbalance  of  active  loads  acting 
on  the  block.  In  accordance  with  the  frictional  hypothesis,  the  maximum 
frictional  resistance  obtainable  is  the  block  weight,  W,  times  the 
coefficient  of  friction,  u . Without  loss  of  generality,  let  W = and 
U = 1 so  that  maximum  frictional  resistance  is  given  bv  the  constant  a . 

y 

The  equilibrium  equation  of  the  block  is  (a  - a ) - n.  *=  0,  where 

p f 

h-i-v 

The  active  load,  a,  shall  be  called  the  "applied  stress"  am!  is 

transmitted  to  the  block  by  the  linear  elastic  spring,  E^,  rep)  i.ting 

the  elastic  portion  of  the  model.  Clearly,  the  applied  stress  is  related 

to  the  elastic  strain,  c , by  o = E c . The  other  active  load,  a , 

e - e e P 

shall  be  termed  the  "plastic  tracking  stress"  and  is  related  to  the 


plastic  strain,  , bv  the  relation  c = E c . 

P p p p 


The  spring  modulus,  , 
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may  be  assumed  constant  or  considered  a function  of  plastic  work  or 
plastic  strain.  Lastly,  it  is  observed  that  the  model  implies  that  the 
total  strain  equals  the  elastic  plus  plastic  strain;  i.e.,  e = z + e . 

6 p 

With  these  ground  rules,  it  is  instructive  to  study  the  model  for 
one  cycle  of  loading  over  the  nominal  range  -2.0 < a < 1.5<J  . Begin- 
ning with  a virgin  material  specimen,  a tensile  load  is  applied.  Within 
the  load  range  0 < a < a , the  model  responds  elastically  as  denoted  by 
the  segment  OA  in  the  stress-strain  diagram  in  Figure  2.  Upon  further 
load  increase  (a  > c^),  the  maximum  frictional  resistance  is  exceeded 

(a r = a ) , and  the  block  moves  an  amount  c such  that  equilibrium  is 
f y P 

maintained;  i.e.,  a = a - a . This  response  implies  plastic  deforma- 

p y 

tion  is  occurring  and  is  denoted  by  the  segment  AB  in  Figure  2.  The 
shape  of  segment  AB  is  dictated  by  the  definition  of  E . For  example, 
if  Ep  is  assumed  zero  (or  nearly  zero),  the  model  would  depict  an  elastic- 
perfectly  plastic  material,  and  AB  would  be  a horizontal  line.  If  E^  is 
assumed  a positive  constant,  AB  becomes  the  upper  portion  of  a bilinear 
stress-strain  curve  with  a slope  equal  to  E E / (E  + E ).  Lastly,  if  E 
is  assumed  to  be  variable,  say  a function  of  plastic  work,  then  a variety 
of  shapes  is  possible  as  suggested  in  Figure  2. 

After  having  subjected  the  model  to  the  maximum  load,  o = o^,  the 
load  is  decreased,  which  produces  the  intriguing  result  that  unloading 
is  linear  elastic;  i.e.,  no  block  movement,  only  elastic  spring  movement. 
If  this  characteristic  of  the  model  is  not  obvious  to  the  reader,  recall 
that  Oj.  is  a passive,  resistance,  and  reverse  movement  of  the  block 
cannot  occur  until  the  frictional  resistance  is  fully  reversed  to  its 
maximum  value,  o^.  Clearly,  elastic  unloading  of  this  model  will  always 
have  a stress  range  of  2a^  as  typified  by  the  segment  BC  in  Figure  2. 

Upon  further  unloading,  the  block  moves,  which  causes  plastic  straining 
in  the  reverse  direction;  this  is  denoted  by  the  segment  CD.  As  before, 
the  shape  of  CD  is  dependent  on  the  definition  of  E^.  Lastly,  the  load 
cycle  is  completed  by  reversing  the  load  to  zero,  which  produces  the 
elastic  response  DE  shown  in  Figure  2. 
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Kinematic  re-entry 
Universal  re-entry  (variable) 

Isotropic  re-entry 


Figure  2.  Stress-strain  path  for  kinematic,  isotropic 
and  universal  surface  hardening. 
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Obviously  this  stress-strain  path  of  the  model  closely  resembles 
actual  experimental  data  of  many  ductile  materials  loaded  into  the 
plastic  range.  Consequently,  it  is  reasonable  to  describe  the  concepts 
of  plasticity  using  the  model  as  a visual  aid. 

A "yield  condition"  is  described  by  the  scalar  function,  F,  which 
predicts  whether  or  not  the  current  stress  state  in  the  plastic  range, 
that  is,  on  the  yield  surface.  With  a little  ingenuity,  such  a function 
can  be  found  for  the  one-dimensional  model  by  recalling  plastic  deforma- 
tion (block  movement)  can  only  occur  when  the  absolute  value  of  the  net 
active  loads  has  mobilized  the  maximum  frictional  resistance,  i.e., 


I i 


3 


i 

I 


F(c,op) 


o 

V 


( 2-4 ) 


Clearly,  F can  never  be  greater  than  zero,  since  this  would  violate  the 
equilibrium  of  the  model.  However,  F can,  and  does,  equal  zero  when  the 
active  loads  have  mobilized  full  frictional  resistance,  and  yielding  is 
occurring  or  about  to  occur.  Lastly,  if  F is  less  than  zero,  the  active 
loads  produce  a frictional  resistance  less  than  maximum,  implying  the 
stress  state  is  in  the  elastic  range.  In  short,  F provides  the  means  of 
determining  whether  or  not  the  current  state  of  stress  is  in  the  plastic 
domain  (on  the  yield  surface)  or  in  the  elastic  domain  (within  the  yield 
surface) . 

With  the  above  insights  in  mind,  it  is  possible  to  generalize  the 
one-dimensional  yield  condition  to  a mult.. . .mens ional  stress  state  yield 
condition  (kinematic)  as  follows: 

F(o,o  ) = f(o  - o ) - k (2-5 ) 

~ ~P  ~P  o 


1 

| 


where 


f(2  ' °p} 

k 

o 


= loading  function 
= yield  parameter  (constant) 


As  before,  F is  the  yield  condition  with  the  following  properties: 


F > 0 Impossible.  A violation  of  equilibrium, 

i.e.,  inadmissible  stress  state.  (2-6a) 

F = 0 Implies  current  stress  state  is  on  yield 
surface,  i.e.,  in  plastic  domain.  (All 
states  of  stress  o satisfying  F = 0 form 

the  current  yield  surface.)  (2-6b) 

F < 0 Implies  current  stress  state  is  within 

yield  surface,  i.e.,  in  elastic  domain.  (2-6c) 


The  loading  function  f(o  - o ) is  the  single  most  important  concept 

- -»p 

in  plasticity  theory.  It  is  a positive-valued  scalar  function  that 
measures  the  magnitude  of  a "select  portion"  of  its  argument  a - a that 

p 

is  responsible  for  plastic  yielding.  In  the  one-dimensional  example 

2 

(o  - Op)  , which  implies 

the  entire  argument  is  responsible  for  plastic  yielding.  However,  in  a 
general  multidimensional  stress  state,  the  vector  argument  of  the  load- 
ing function  allows  a great  deal  of  freedom  in  choosing  the  form  of  the 
loading  function.  For  example,  some  materials,  such  as  ductile  metals, 
exhibit  plastic  responses  due  to  shearing  stresses,  but  show  no  plastic 
response  when  subjected  to  hydrostatic  stresses.  Thus,  it  would  be 
appropriate  to  define  a loading  function  that  increases  with  shear 
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stress  but  not  with  hydrostatic  stress,  such  as  the  second  invariant  of 
deviatoric  stress  2^.  Other  materials,  such  as  soils,  demonstrate 
plastic  responses  due  to  various  combinations  of  the  stress  invariants. 

The  literature  is  full  of  a wide  variety  of  proposed  loading  func- 
tions for  different  materials.  Several  common  loading  functions  are 
presented  in  Appendix  A.  For  now,  it  is  only  important  to  appreciate 
that  the  loading  function  acts  like  a filter  in  the  sense  it  only 
increases  or  decreases  its  value  for  some  predefined  "select  portion"  of 
the  stress  vector  considered  to  be  responsible  for  plastic  yielding. 

All  other  components  of  stress  are  filtered  out  and  do  not  influence  the 
value  of  the  loading  function. 

To  better  appreciate  the  significance  of  the  loading  function, 
consider  Figure  3.  Assume  F = 0,  which  implies  the  current  state  of 
stress  is  on  the  yield  surface.  Also  suppose  an  arbitrary  stress 
increment,  A o,  is  applied  to  the  current  stress  state.  V.'ill  this  cause 
plastic  deformation?  To  answer  this,  it  is  merely  necessary  to  know  if 
the  stress  increment  has  any  components  in  the  direction  of  the  outward 
normal  of  the  loading  function,  i.e.,  in  the  direction  that  the  yield 
surface  moves.  The  dot  product  of  the  stress  increment  with  the  loading 
function  gradient  (which  is  in  direction  of  outward  normal)  provides  a 
simple  test  with  three  possible  results: 


(1 


Ao  > 0 


Plastic 
portion 
surf ace 


loading;  i.c.,  at  least  some 
of  Ac  is  col  inear  with  outward 
normal.  Surface  moves  out. 


C2-7a) 


0 Neutral  loading,  no  plastic  deformation 
occurs;  i.e.,  Ac  is  in  the  tangent  plane 
of  the  loading  surface.  Surface  does  not 
move.  (2- 7b) 
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f • Ao  < 0 

~o  ~ 


Elastic  unloading  occurs;  i.e.,  at  least 
some  portion  of  do  is  colinear  with  the 
Inward  surface  normal.  Surface  does  not 
move . 


(2-7c) 


where  f = 3f/3o  is  the  loading  function  gradient,  and  is  in  the  direc- 
tion of  the  outward  normal. 

Clearly,  the  importance  of  the  loading  function  cannot  be  over- 
stressed since  it  dictates  what  stress  increments  promote  plastic  yield- 
ing. 

The  yield  parameter,  , in  Equation  2-5  is  a material  dependent 
constant  and  may  be  geometrically  interpreted  as  the  "radius"  of  the 
yield  surface.  The  value  of  k^  may  be  determined  from  any  type  of 
standard  laboratory  test  by  evaluating  the  loading  function  at  the 
stress  state  producing  initial  yielding,  i.e.,  kQ  = f(n^),  where  is 
the  stress  vector  producing  initial  yielding. 

In  general,  "surface  hardening"  is  a simple  concept,  as  it  merely 
pertains  to  keeping  track  of  the  "radius"  and  "center"  of  the  yield 
surface  (not  to  be  confused  with  hardening  rule  for  H'  discussed  later). 
For  the  kinematic  model  under  discussion,  the  radius  of  the  yield  sur- 
face remains  constant,  i.e.,  o for  a one-dimensional  model,  and  k for 

y o 

the  general  model.  However,  the  "center"  translates  every  time  plastic 
deformation  occurs.  The  amount  of  translation  can  be  deauced  by  con- 
sidering the  total  derivative  of  the  yield  condition. 

To  see  this,  first  consider  the  one-dimensional  model.  Initially, 

= o - u,  and  the  radius  is  ' with  the  center  at  the  origin.  As  soon 

i sf led  , 

i.e.,  F = 0 and  the  total  derivative,  dF  = 0,  must  also  be  satisfied. 

That  is,  from  Equation  2-4: 


0,  and  the  radius  is  ' with  the  center  at  the  origin. 

P y 

as  the  load  o reaches  the  yield  stress,  the  yield  condition  is 


= 0 = 


(2-bl 
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In  the  above  Ao  and  Ao  are  actuallv  differentials  do  and  da  ; however, 

P P 

throughout  this  writing  the  incremental  symbol  "A"  will  be  used  to 

emphasize  numerical  approximations.  Equation  2-8  simply  states  that  any 
increment  Ao  of  the  applied  stress  is  immediately  followed  by  an  equal 
increment  Aa^  of  the  plastic  tracking  stress.  This  implies  the  "center" 
of  the  elastic  range  shifts  an  amount  Ac^ , while  the  "radius"  of  the 
elastic  range  remains  a constant  . Noting  is  the  sum  of  its 
increments,  = lAa^,  a geometric  interpretation  of  the  plastic  track- 
ing stress  can  be  offered.  Namely,  a is  the  current  center  of  the 
elastic  range  as  it  translates  up  or  down  the  stress  axis.  This  concept 
is  demonstrated  in  Figure  4,  wherein  the  center  and  radius  of  the  elastic 
range  are  shown  at  various  load  points  corresponding  to  Figure  2. 

This  geometric  interpretation  of  the  plastic  tracking  stress  also 
applies  to  multiaxial  stress  states  in  that  the  yield  surface  translates 
rigidly  (constant  radius  kQ)  in  stress  space  with  its  "center"  located 
by  0^  as  suggested  in  Figure  5. 

To  prove  this  for  the  general  kinematic  case,  it  is  necessary  to 
assume  the  direction  of  Ao^  during  yielding.  In  the  one-dimensional 
case,  this  was  not  necessary  since  only  one  direction  was  possible. 

Prager  [18]  made  the  assumption  that  Ao^  is  in  the  direction  of  the 
outward  normal  of  the  loading  function,  i.e.,  in  the  direction  of  its 
grad ient . 

In  light  of  Equation  2-7b,  Prager's  assumption  has  great  intuitive 

appeal  since  it  implies  Aa  is  in  the  direction  cf  the  select  portion  of 

~P 

the  applied  stress  that  causes  plastic  yielding.  Therefore,  it  is 
assumed : 


Ao  = C f 

~P 


(2-9) 


where  f 

“V  ft 


Sf/.Tj,  loading  function  gradient  with  respect  to 


= scalar  constant  to  be  determined 
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Track  of  Elastic  Stress  Range,  ola 


- Oy*/Oy  at  E 


Figure  4.  One-dimensional  kinematic 
surface  hardening. 


Figure  5.  Multidimensional  kinematic  surface  hardening. 


The  value  of  C Is  obtained  by  taking  the  total  derivative  of  the 
yield  condition  (Equation  2-5)  which  must  be  zero  during  yielding,  i.e.. 


dF  = 0 = f 

-.o 


Ao  + f 

~ ~o 


(2-10) 


where  fG  = 3f/8o  , loading  function  gradient  vector  with  respect  to  o . 
~ p ~P  ~P 

Because  the  argument  of  the  loading  function  has  the  special  form, 

o-o  , it  is  trivial  to  show  that  f = -f„  for  anv  function  f.  There- 

~ ~P  ~~  ~°o 

fore,  Equation  2-9  can  be  combined  with  2-10  to  give: 


C = (fj  * Ao)/(f  T • f ) 


Returning  to  Equation  2-9,  the  desired  relationship  is  obtained: 


, ,aT  .a 

Ao  = (n  Ao)n 

~p  


(2-11) 


n = f /J  f T f 


(2-12) 


where  n is  the  unit  outward  normal  of  the  loading  function. 

Equation  2-11  provides  the  rule  for  determining  the  increment  of 

plastic  tracking  stress  for  every  applied  stress  increment  Ac.  The 

total  sum  of  the  plastic  tracking  stress  increments,  i.e.,  o = EAo  , 

-P  ~P 

provides  the  coordinates  for  tracking  the  "center"  of  the  yield  surface 
as  was  suggested  in  Figure  5. 
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This  completes  the  general  development  for  the  kinematic  yield 
condition  and  surface  tracking  rules.  Next,  isotropic  hardening  is 
examined . 


Isotropic  Model 

Figure  lb  is  a one-dimensional  representation  of  an  elastic-plastic 

model  with  isotropic  hardening.  All  of  the  ground  rules  previously 

described  for  the  kinematic  model  still  hold.  The  only  difference  is 

that  the  spring  has  been  moved  to  act  vertically  on  the  block  as 

shown.  Furthermore,  it  is  imagined  some  servo-mechanism  compresses  this 

spring  an  amount  je  | everytime  the  block  moves  a distance  in  either 

direction.  The  accumulation  of  all  plastic  strain  movements  is  denoted 

by  e * ; i.e.,  c * = lie  |. 

P P ' P 

The  significant  consequence  of  this  new  model  is  that  the  load  c 

P 

has  been  transformed  from  the  status  of  an  active  load  into  a passive 

frictional  resistance  and  is  hereby  renamed  as  o^*.  Accordingly,  the 

maximum  frictional  resistance  is  no  longer  constant,  but  rather  is  given 

bv  Of  =o  + a *,  where  o * = E e *.  Equilibrium  of  the  model  is 
1 max  >’  P P P P 

given  by  o - cf  = 0,  where  o < Of 

1 * max 

With  the  above  in  mind,  it  is  instructive  to  trace  the  performance 

of  the  new  model  through  a loading  cycle  as  was  previously  done  for  the 

kinematic  model.  Beginning  with  a virgin  specimen  and  referring  back  to 

Figure  2,  the  initial  elastic  loading  follows  the  same  path  0A  as  the 

kinematic  model.  Furthermore,  the  identical  curve  AB  is  traced  with  a 

loading  above  o . This  is  because  c * = c , which  implies  o * = o , 
y P P P P 

and,  hence,  both  models  have  identical  equilibrium  equations.  At  poii  t 

B unloading  begins.  The  block  cannot  move  backward  until  the  applied 

stress  o completely  reverses  the  current  maximum  frictional  resistance 

(°f  = ou).  Therefore,  elastic  unloading  continues  to  point  C' , where 

max  b 

c = ~Jg-  With  additional  negative  loading  the  curve  C'D'  is  traced. 

Note  C'D’  has  the  same  shape  as  CD;  however,  it  is  rigidly  shifted  along 
the  elastic  unloading  path  an  amount  C'C  = 2(o  - c ). 

D V 
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It  is  evident  that  the  kinematic  and  isotropic  models  will  produce 

extraordinarily  different  results  after  a few  loading  cycles  even  though 

E and  E are  the  same  for  both  models.  However,  the  models  will  give 
P e 

identical  results  as  long  as  the  load  is  monotonically  increasing. 

The  yield  condition  for  the  one-dimensional  isotropic  model  can  be 
written  as: 


p (o,0v*) 


(2-13) 


where  o * - 3 + o *,  the  vield  stress  parameter. 

y >’  P 

Generalizing  the  isotropic  yield  condition  to  multidimensional 
stress  states  gives: 


F(o,k*)  » f(o)  - k* 


( 2-  ] A ) 


where  f(,a) 
k* 


loading  function 

yield  parameter  (increases  with  plastic  def evtnation) 


As  before,  the  yield  condition,  F,  denotes  whether  the  current 
state  of  stress  is  in  the  elastic  or  plastic  domain,  and  has  the  proper- 
ties given  by  Equations  2-6a,  2-6b.  and  2-6c. 

The  concept  of  the  loading  function  is  the  same  as  described  for 
the  kinematic  model,  and  it  has  the  properties  given  by  Equations  2-7a, 
2-7h,  and  2-7c.  The  only  difference  is  that  the  argument  is  composed  of 
the  stress  vector  o rather  than  the  stress  vector  difference  o - o ; 
consequently,  the  "center"  of  the  yield  surface  remains  fixed  at  the. 
origin  of  stress  space. 


1 
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The  "radius"  of  the  yield  surface  is  measured  by  the  yield  param- 
eter k*  which  increases  with  plastic  deformation.  A pictorial  represen- 
tation of  isotropic  surface  hardening  for  the  one-dimensional  load  cycle 
is  shown  in  Figure  6,  whereas  Figure  7 illustrates  isotropic  hardening 
in  a multidimensional  stress  state. 

To  determine  the  amount  the  yield  surface  expands  for  every  applied 
stress  increment  causing  plastic  deformation,  the  total  derivative  of 
Equation  2-14  is  set  to  zero  to  give: 


dF  = 0 = f T ■ A a - Ak*  (2-13) 

Equation  2-15  provides  the  surface  hardening  rule  for  updating  the 

"radius"  of  the  yield  surface  after  every  increment  of  applied  stress. 

Specifically,  k*  = k + F.Ak*,  where  k is  the  initial  vield  parameter 

o o 

(radius)  and  Ak*  = f^T  • Ao. 

Universal  Model 

The  kinematic  and  isotropic  responses,  as  displayed  in  Figure  2, 
represent  two  extreme  predictions  for  re-entry  into  the  plastic  domain. 
Most  ductile  materials  exhibit  a re-entry  point  "C"  somewhere  between 
the  points  C and  C'.  Clearly,  it  would  be  desirable  to  choose  a surface 
hardening  model  that  allows  some  flexibility  in  selecting  the  ’’radius" 
and  "center"  of  the  yield  surface. 

The  "universal  surface  hardening"  concept  provides  this  flexibility 
by  taking  a linear  combination  of  the  kinematic  and  isotropic  m.v'els. 
Figure  lc.  illustrates  this  concept  with  a one- di mens iona  1 model,  wherein 
3 is  a weighting  parameter  (0  < 8 < 1)  applied  to  the  vertical  load, 

0^*,  and  the  remaining  weight  1-6  is  applied  to  the  horizontal  load,  . 
If  6=0,  the  universal  model  reduces  to  the  isotropic  form;  if  6 = 1, 
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Load  Points  from  Figure  2. 


Figure  6.  One-dimensional  isotropic  surface  hardening. 


Fieure  7.  Multidimensional  isotropic  surtace  hardening 


the  kinematic  form  is  attained.  Any  other  value  of  £ combines  the 
i ' characteristics  of  both  models  so  that  both  radius  and  center  of  the 

j yield  surface  change  during  plastic  deformation.  For  example,  the 

| universal  model  will  re-enter  the  plastic  range  at  point  C"  in  Figure  2 

| :■  is  chosen  as  the  ratio:  £ = CC"/CC' . 

| The  ground  rules  for  this  one-dimensional  model  are  the  same  as 

i stated  for  the  kinematic  and  isotropic  models;  therefore,  during  initial 

monotonic  loading,  the  stress-strain  curve  of  this  model  follows  the 
same  path  OAB  in  Figure  2 for  any  value  of  £.  However,  after  elastic 
! unloading  and  re-entry  into  the  plastic  domain  at  point  C"  the  stress- 

strain  path  for  all  future  load  cycles  is  a unique  function  of  £. 

(Note,  for  added  generality,  B could  be  specified  as  a function  of  the 
! number  of  load  cycles,  6 = B(n),  thereby  allowing  the  mode]  to  shift 

: from  isotropic  to  kinematic  or  vice-versa  during  load  cycling.) 

Equilibrium  of  the  universal  model  can  easily  be  deduced  as: 


(a  - B a ) - a = 0 

P * 


where 


a + (1  - S)o  * 
y P 


o 


f 


max 


I 


The  corresponding  one-dimensional  yield  condition  is  given  by: 


F(-,o  .3  *)  = V (o  - 5 r - 5 * 

p y » p y 


(2-lb) 
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where 


o = 6 a 

P P 

o * = o + (1  - 6)o  * 

y y p 


Generalizing  the  universal  yield  condition  to  multidimensional 
stress  states  gives: 


F (o,o  ,k*)  = f(o  - o ) - k* 

~ ~P  ~ ~P 


(2-17) 


where  = 6 o^,  "weighted"  plastic,  tracking  stress 

k*  = k + (1  - S)k*,  weighted  yield  parameter 


As  always,  F is  the  yield  condition  with  the  properties  denoted  in 
liquations  2-6a,  2-6b,  and  2-6c,  and  f is  the  specified  loading  function 
whose  gradient  can  be  used  for  the  plastic  yielding  test  denoted  by 
Equations  2-la,  2-7b,  and  2-7c. 

The  vector  5^  and  the  scalar  ky<  specify  the  origin  anu  the  radius, 

respectively,  of  the  universal  yield  surface.  Figure  8 demonstrates  the 

concept  of  universal  hardening  for  the  one-dimensional  load  cycling, 

while  Figure  9 portrays  an  example  evolution  of  the  yield  surface  in 

multidimensional  stress  space. 

To  determine  the  amount  of  radius  expansion  Ak*  and  center  shift 

An  for  any  applied  stress  increment  causing  plastic  deformation  (F  = 0) , 

the  total  derivative  of  Equation  2-17  could  be  set  to  zero  to  permit 

- T 

solving  for  do  and  dk*  in  terms  of  f and  do.  However,  it  is  tar 

P *m\J  *** 

simpler  to  directly  vise  the  definitions  of  c and  k*  from  Equation  2-17 
to  determine  the  derivatives;  i.e.,  Ak*  = (1  - B)Ak*  and  An  = 6 An  . 

~p  ~p 

With  the  aid  of  Equations  2-11,  2-12,  and  2-15,  the  increments  of  the 
radius  and  center  of  the  universal  vield  surface  are: 


Thus,  the  universal  surface  hardening  rule  may  be  written  as: 


( 

i 


(2-20) 


(2-21) 


where  k*  is  the  current  "radius"  of  the  vield  surface,  and  o is  the 

~P 

current  vector  from  the  origin  of  stress  space  to  the  "center"  of  the 
yield  surface. 

Flow  Rule 

The  purpose  of  a flow  rule  is  to  define  increments  of  plastic 
strain  accompanying  an  increment  of  applied  stress  during  plastfe 
yielding  (F  = 0).  For  any  and  all  of  the  one-dimensional  models  this  is 
a straight-forward  application  of  the  ground  rules,  and  it  is  easy  to 
deduce  the  flow  rule  to  be: 


Figure  8.  One-dimensional  universal  hardening  with 
6=1/2. 


where  E ' is  the  tangent  value  of  the  uniaxial  curve  relating  total 
plastic  strain  to  total  stress  as  illustrated  in  Figure  10. 

When  the  flow  rule  is  generalized  to  multidimensional  stress-strain 
states,  the  development  is  not  quite  as  simple.  The  complication  is  due 
to  the  fact  that  now  the  plastic  strain  increment  is  a vector;  there- 
fore, it  must  be  defined  with  respect  to  direction  as  well  as  to  magni- 
tude. (Note,  in  the  one-dimensional  idealization,  direction  and  mag- 
nitude were  one.) 

Thus,  for  the  general  case,  the  plastic  strain  increment  vector  may 
be  described  by 


a£ 

-P 


dX  m 


(2-23) 


where  m is  some  directional  vector  (in  stress  or  strain  space)  of  unit 
magnitude,  and  dX  Is  a scalar  denoting  the  magnitude  of  the  plastic 
strain  increment.  To  determine  dX  and  m require  assumptions  similar  to 
the  assumptions  made  for  the  plastic  tracking  stress  vector. 

The  first  assumption  is  based  on  the  observation  that  the  magnitude 
of  the  plastic  strain  increment  should  be  proportional  to  the  magnitude 
of  the  applied  stress  increment  which  is  colinear  with  the  outward 

A 

normal  (n)  of  the  loading  function;  i.e.. 


,,  1 ,aT  , 

d ' = H~'  ~ j 


(2-24) 


The  inner  product,  A o,  is  the  scalar  measure  of  stress  that  pushes 

the  yield  surface  outward,  and  H'  is  a hardening  coefficient  (modulus) 
determined  from  experimental  tests.  H*  will  be  discussed  further  in  the 
hardening  rule. 
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Figure  10.  Shape  of  stress-plastic  strain  curve. 
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With  regard  to  the  direction  m of  the  plastic  strain  increment,  two 
assumptions  are  available.  By  far  the  most  common  is  the  so-called 
"associative  flow  rule  assumption."  Here,  the  idea  is  that  plastic 
straining  occurs  in  the  same  direction  as  the  outward  normal  of  the 
loading  function.  Although  this  may  be  aesthetically  pleasing,  it  is 
still  only  an  assumption.  Other  plastic  potential  functions  (loading 
functions)  could  just  as  well  be  used  to  predict  the  direction  of  the 
j.  lastic  strain.  In  such  cases,  the  procedure  is  termed  "nonassoc iative 
flcv>  rule  assumption." 

In  short,  the  direction  m has  two  possibilities: 

m * n associative  flow  law  (2~25a) 

m j*  n nonassociative  flow  law  C 2—  25b) 


More  on  nonassociative  flow  rules  is  given  in  Appendix  B. 

Inserting  Equation  2-24  into  Equation  2-23  and  with  the  understand- 
ing of  Equation  2-25,  the  general  flow  rule  is: 


At  = (nT  Ao)m  (2-26) 

~p  ri  '■w  -v 

Hardening  Rule 

The  hardening  coefficient,  H',  is  the  only  quantity  that  remains  to 
be  defined.  The  value  of  H'  must  be  determined  from  experimental  results 
consistent  with  Equation  2-26.  That  is,  H'  is  dependent  on  n and  m as 
well  as  c and  c . As  an  example,  suppose  it  is  decided  to  determine  H' 


41 


by  means  of  a simple  tension  test.  Let  the  curve  in  Figure  10  represent 
axial  stress  versus  plastic  strain  (i.e.,  elastic  strains  have  been 
subtracted  from  total  strains).  Since  all  components  of  stress  are  zero 
except  cr^,  the  inner  product  in  Equation  2-26  is  simply  Ao^n^.  Also, 
Equation  2-26  is  valid  for  every  component  of  plastic  strain,  in  parti- 
cular, it  is  valid  for  Atp^.  Therefore,  H’  = m^n^ (Ac Atp^ ) , or  using 
Equation  2-22  the  hardening  coefficient  determined  from  a tensile  (or 
compression)  test  is: 


H’ 


m. 


E ' 
P 


(2-27) 


For  the  special  case  where  a Von  Mises  loading  function  *~d  an  associative 

flow  rule  is  assumed,  n,  = m.  = 2/\  6,  and  H'  = 2/3  E '. 

I i » p 

Furthermore,  it  is  emphasized  that  the  determination  of  H'  is  not 
restricted  to  tensile  tests.  Any  type  of  experimental  data  can  be  used 
to  determine  H'  providing  Equation  2-26  is  used  consistently. 

In  the  event  H'  is  constant  throughout  the  loading  path,  no  addi- 
tional assumptions  are  necessary,  and  the  flow  rule  is  well  defined  for 
all  stress  states.  However,  in  general,  H'  varies  throughout  the  load 
path  as  implied  in  Figure  10.  This  gives  rise  to  the  need  of  a "harden- 
ing rule,"  wherein  H'  is  assumed  to  be  a function  of  some  plastic  response. 
The  first  inclination  is  to  consider  H'  a function  of  the  magnitude  of 
plastic  strain,  i.e.,  H'  = r.  °,  where  t:  L is  plastic  strain  magnitude 

o r~\ r p p 1 

given  by  e ^ ^ • < . Although  this  assumption  is  adequate  tor  a 

monotonic  radial  loading,  for  general  loadings  it  produces  results 

unrepresentative  of  actual  material  behavior.  For  example,  consider  the 

response  of  the  one-dimensional  kinematic  model  shown  in  Figure  2.  If 

E ' = E '(c  °) , then,  during  the  reverse  loading  phase,  the  value  of  t ( 

P P P P 

reduces.  This  results  in  an  increased  value  of  E ' , and  the  shape  of 

P 

segment  CD  is  concave  up  rather  than  concave  down  as  desired. 
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To  avoid  this  anomaly,  the  hardening  rule  must  be  based  on  some 
nondecreasing  measure  of  plastic  response.  Two  such  measures  are 
commonly  employed:  plastic  work  and  "total"  plastic  strain.  The  former 
Is  known  as  the  "work  hardening  rule"  and  simply  asserts  that  H'  is  a 
function  of  the  plastic  work  per  unit  volume,  W^: 


H'  = H!(W  ) 
P 


(2-28) 


where  W 

P 


dc 

~P 


The  second 
assumed  that  H' 
magnitude,  £p*: 


method  is  known  as  "strain  hardening,"  wherein  it  is 
is  a function  of  the  "total  path"  of  the  plastic  strain 


H'  = H'(c  *) 
P 


(2-29) 


where  t * = ^ Ac  ° 

P t-t  P 

Note:  e * is  not.  the  same  measure  as  c °,  since  c * increases  with 

P o P P 

every  increment  Ac^;  however,  is  the  measure  of  current  plastic 

strain  magnitude  regardless  of  the  path. 

The  choice  of  one  hardening  rule  over  the  other  is  largely  t matter 
of  computational  convenience.  By  either  rule,  H'  is  generally  not  known 
as  a continuous  function,  but  rather  as  discrete  points  of  the  plastic 
work  or  total  plastic  strain  as  suggested  in  Figure  11.  It  is  a trivial 
matter  to  store  the  discrete  points  in  the  computer  and  interpolate  to 
find  the  current  value  of  H'. 


INCREMENTAL  PLASTICITY  RELATIONS 

The  four  basic  concepts  of  plasticity  have  been  presented.  The 

remaining  job  is  to  pull  the  pieces  together  to  find  the  matrix  quantity 

D denoted  in  Equation  2-1. 

~ep  n 

Beginning  with  the  combination  of  Equations  2-2  and  2-3,  Aa  is 
given  as: 


Ao  = D (Ac  - Ac  ) 

"n*  om  G ~P 


(2-30) 


Taking  the  dot  product  of  Equation  2-30  with  respect  to  the  unit  gra- 

at 

dient  of  the  loading  function  u , and  replacing  Ac  by  Equation  2-26, 

aT  ~ 

the  scalar  quantity  (n  Ao)  can  be  determined  as: 


AT  , 

n A o 


nT  D Ac 
- «e 

1 + (nT  D m)/H' 
~ -e  ~ 


(2-31) 


Using  Equation  2-31  in  Equation  2-26,  the  plastic  strain  increment 
is  related  to  the  total  strain  increment  by: 


C D Ac 

-p  me  ~ 


(2-32) 


where 


H'  + (n^  D m) 
*»  ™e  ~ 
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LastLy,  inserting  Equation  2-32  into  Equation  2-30,  the  desired  stress- 
strain  relationship  is  achieved: 


(D-D  ) Ae 
*e  «.p 


(2-33) 


where  D 

-P 


D C D 
■»e  «*p  «*e 


Equation  2-33  is  the  key  incremental  stress-strain  relationship  for  use 
in  elastic-plastic  boundary  value  problems.  The  plasticity  matrix,  D^, 
is  dependent  on  the  variables  H'  and  n (and  m for  nonassociative  law). 
Traditionally,  these  quantities  are  updated  at  the  beginning  of  each 
load  step  and  are  assumed  to  be  constant  over  the  load  increment.  These 
concepts  will  be  discussed  later  when  the  plasticity  relationships  are 
incorporated  into  a finite  element  formulation.  For  now,  the  plasticity 
constitutive  theory  is  reviewed. 


SUMMARY  OF  CONSTITUTIVE  THEORY  OF  PLASTICITY 


For  easy  reference,  the  pertinent  plasticity  relationships  for 
universal  hardening  are  listed  below. 

The  universal  yield  condition  is: 


F(c,5  , k*)  = f(o  -a  ) - k* 

~ ~P  ~ ~P 


(2-34) 


F is  the  yield  function,  such  that  F = 0 implies  yielding,  c is  the 

total  stress  vector,  a is  the  plastic  tracking  stress  locating  the 

~P 


*'*■ ■Lij.i'j ^fttiiliyuf  a iiiMi 


"center"  of  the  yield  surface,  and  k*  is  a scalar  stress  measure  defining 
"radius"  of  the  yield  surface.  The  function  f is  the  prescribed  loading 

function. 

The  surface  hardening  law  is  given  by: 


Ak*  = (1  - 6) f T Ao  (2-35) 

** 

Aa  = B(m  • Ao)n  (2-36) 

~ ~ ~ 


where 


(2-37) 


(2-38) 


Equations  2-35  and  2-36  give  the  rule  for  keeping  track  of  the  yield 
surface,  i.e.,  k*  = kQ  + EAk*  and  5p  = EA£p.  , is  the  universal  harden- 
ing parameter  (0  < g < 1),  where  6=0  implies  isotropic  hardening, 

6 = 1 implies  kinematic  hardening,  and  other  values  of  6 imply  a linear 
combination  of  kinematic  and  isotropic  hardening.  In  Equation  2-37  f(j 
is  the  gradient  of  the  loading  function  and  is,  by  definition,  in  the 

A 

direction  of  the  outward  normal  of  the  yield  surface.  Accordingly,  n in 
Equation  2-38  is  the  unit  outward  normal. 

The  hardening  rule  is: 


H'(c  *) 
P 


(2-39) 
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(2-40) 


c * - Ac  * 

F P 

Ac  * %/Ac  T . Ac  (2-41) 

P \ ~P  ~P 


Here  H'  is  known  as  the  hardening  coefficient  to  be  used  in  the 
flow  rule.  It  is  assumed  H'  is  a function  of  the  accumulated  plastic 
strain  norm,  c^*.  The  function  H'(c  *}  is  prescribed  material  data. 
The  flow  rule  is  given  by: 


AC 

-P 


1 

H 


Cn^  Ao ) m 


(2-42) 


The  direction  vector,  m,  is  the  direction  cf  plastic  straining. 

If  in  “ n,  Equation  2-42  is  termed  an  "associative"  flow  rule.  If  m is 
defined  from  some  other  potential  function.  Equation  2-42  is  termed 
"nonassoc iative. " In  either  case,  the  flow  rule  can  be  used  to  deter- 
mine the  incremental  plastic  constitutive  law: 


Aa  = (D  - D ) Ac 

■»e  H p ~ 


(2-42) 


D 

«P 


rv  A aT  ,, 

D m n 'J 

m C ~ — m€ 

H'  + (nT  D m) 

-v  >»e  - 


(2-44) 


Equations  2-34  through  2-44  summarize  the  incremental  laws  of  plasticity. 
In  the  next  section  these  laws  will  be  incorporated  into  a finite  element 
formulation. 
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PLASTICITY  FINITE  ELEMENT  FORMULATION 


From  the  previous  finite  element  development,  the  general  equilib 
rium  equation  (Equation  1-6)  was  given  as: 


Aa  dv  = AP 


(2-45) 


In  the  above,  B,  is  the  strain-to-nodal  displacement  matrix,  AP  is  the 
external  load  increment,  and  the  symbol  £ implies  the  ordered  summation 
of  the  element  volume  integrations. 

Inserting  the  constitutive  law.  Equation  2-43,  into  the  equilibrium 
equations,  and  using  Ae  = B Au,  Equation  2-45  can  be  written  as: 


(K  - K )Au  = AP  (2-46) 


K 

-e 


Eft 


D 

we 


B dv 


(2-47) 


K 


Ef 


BT  D B dv 

X ,p  x 


(2-48) 


In  the  above',  Kp  is  the  global  elastic  stiffness  matrix  which  remains 

constant  throughout  the  loading  schedule.  K is  the  global  plastic. 

*P 

stiffness  (reduction)  matrix  which  is  dependent  on  stress  state.  Note, 

K - 0 whenever  the  stress  state  is  within  the  vield  surface.  As  dis- 
ap 

played  in  Equation  2-46,  the  equilibrium  equations  are  in  the  proper 
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form  for  a direct  solution  by  the  tangent  stiffness  method.  This  method 
requires  the  combined  stiffness  matrix  Kg  - to  be  triangular ized  at 
least  once  for  each  load  step  in  which  plastic  deformation  occurs. 
Although  this  method  may  be  inefficient,  many  researchers  prefer  this  to 
the  uncertainties  of  iterative  techniques. 

Alternatively,  to  obtain  the  proper  form  for  the  initial  strain 
method.  Equation  2-46  is  written  as: 


K Au  = AP  + AF 
«e  ~ ~ **p 


(2-49) 


where 


AF  = p / BT  D B dv 
~P  •/  m Mp  ** 


(2-50) 


In  Equation  2-49,  K needs  only  to  be  triangularized  once.  All  non- 

*e 

i inearities  are  introduced  through  the  unknown  force  vector  AF  , and  an 

~P 

iterative  method  is  employed. 

Both  solution  procedures  are  outlined  in  the  following  pages. 


Tangent  Stiffness  Method 


It  is  assumed  the  following  quantities  are  known  at  load  step  "n" 
u , a , and  e , along  with  the  current  plastic  measures  D , k *,  tD  *, 


~n  ~n 


n 


n 


and  nD  . The  objective  is  to  find  the  increments  of  the  above  quantities 
~ 1 n A 

from  load  step  n to  n + 1;  i.e.,  Au,  Ao,  Ac,  etc.,  and  to  update  the 
plastic  measures. 

1.  Determine  load  increment  AP. 


2.  Assemble  K - K from  Equations  2-47  and  2-48,  where  K is 

a>e  »p  *p 

based  on  at  the  beginning  of  the  step. 
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AP,  for  au. 


Solve  (K.  - K ) An  - 

■ e «p  ~ 

Compute  stress  and  strain  increments  for  each  element  (or  each 
integration  point);  i.e., 

a.  At  = B A u 

b.  Ac  = (D  - D ) Ae 

..  »e  «p  ~ 

Monitor  each  element  (or  integration  point)  to  determine  if  the 
current  stress  state  is  in  the  plastic  region  or  elastic  region: 


a.  Compute  F 

- f(o  ,,  - c ) - k * 
~n+l  p n 

b.  If  F < 0, 

stress  state  is  elastic 

. Set  D = 0 

-P 

and  go  to 

Step  8. 

C.  If  F > 0, 

stress  state  is  plastic 

. Go  to  Step 

6. 

Update  yield 

surface  measures: 

a.  o = 

~pn+l 

o + An  , where  Aci  - 

~Pn  ~P  ~P 

« fnT  Ac)n 

b‘  Vfl*  = 

f (o  . - o ) 

-n+1  ~Pn+l 

Update  plastic  flow  law: 

^ f 

a • n - t / 
- 

f ^ where  f is 

evaluated  at 

stress  state 

n + 1 . 

b.  m = n for  associative  flow  rule;  otherwise  compute  m from 
nonassoc i a t ive  function. 

c.  H’  = H’(r.  *),  where  t_  * = * + 


■■  I 


1 


d.  D = D ra  n D / [H*  + (n  D ffi)  ] 

»p  „e  ~ - »e  ~ »e  ~ 

3.  Print  out  desired  results,  and  return  to  Step  1 for  next  load 
increment . 


The  preceding  algorithm  is  a relatively  straightforward  procedure; 
however,  it  is  not  computationally  efficient.  Worse  still,  the  transi- 
tion error  may  require  additional  matrix  triangularizations  within  the 

load  step.  It  is  assumed  the  matrix  D remains  constant  (chord  value) 

*P 

during  a given  load  step.  This  assumption  is  reasonable  providing  the 
element  is  in  the  plastic  region  at  the  beginning  and  end  of  the  load 
step,  or  it  is  in  the  elastic  region  at  the  beginning  and  end  of  the 
load  step.  However,  when  an  element  shifts  from  the  elastic  zone  to 
the  plastic  zone  or  vice  versa  (the  transition  region),  it  must  be 


treated  with  special  care.  In  these  cases,  actually  abruptly  changes 
from  zero  to  some  finite  value  as  the  element  transcends  from  an  elastic 


state  to  a plastic  state,  or  conversely,  from  a finite  value  to  zero  for 
the  reverse  transition. 

To  account  for  these  transitions  in  the  above  algorithm,  two  basic 
approaches  are  available.  First,  and  simplest,  is  to  divide  the  load 
into  sufficiently  small  increments  so  that  the  transition  error  can  be 
ignored.  By  this  method,  the  algorithm  lags  the  transition  responses  by 
one  load  step. 

The  alternative  approach  is  to  modify  is  some  way  to  account  for 
tran>ition  and  resolve  the  problem.  In  the  case  of  the  transition  from 

plastic  to  elastic  zone,  D should  be  set  to  zero  because  all.  unloading 

" P 

is  elastic.  However,  when  the  transition  is  from  the  elastic  to  the 
plastic  zone,  the  first  part  of  the  load  step  is  in  the  elastic  zone  as 
the  stress  path  moves  from  some  point  within  the  yield  surface  to  the 
yield  surface.  The  second  part  of  the  load  step  is  in  the  plastic  zone 
as  the  stress  path  moves  with  the  yield  surface.  To  reflect  the  abrupt 

change  of  D during  the  load  step,  it  is  convenient  to  take  a weighted 

"P 

average  based  on  the  proportion  of  the  load  increment  in  tin  plastic 
zone  as  measured  bv  the  ratio: 


r 


£'2i+l  - T^1  ' V 

l(2i+i  ' 4.^  ”<£1  ‘ 2Pl> 


(2-51) 


Thus,  in  the  algorithm,  each  element  is  checked  to  see  if  it  underwent  a 
transition  phase.  If  so,  the  load  step  is  repeated,  wherein  for  each 
transition  element  is  re-estimated  as  follows,  For  an  unload  transi- 
tion, set  D = 0.  For  loading,  the  transition  ratio,  r,  is  computed 

• p 

from  Equation  2-51,  and  the  stress  state  at  the  yield  surface  is 
o = + (1  - r )Sa.  Next,  D is  calculated  at  the  stress  state  o and 

a»  -%,  * 1 ^ tm  p *•» 

then  reduced  hv  the  factor  r.  The  reduced  value  of  D is  used  for 

»P 

calculating  the  plastic  stiffness  K , and  the  algorithm  proceeds  as 

»P 

be  fore . 


The  computational  inefficiencies  of  the  tangent  stiffness  method 
are  clearly  evident,  since  the  combined  stiffness  matrix  must  be  assem- 
bled and  triangular ized , not  only  for  each  load  step,  but  also  again 
within  a load  step  when  element  transitions  occur.  The  advantages  of 
the  initial  strain  method  are  demonstrated  in  the  next  section. 


Initial  Strain  Method 


Equation  2-49  represents  the  governing  equations  for  the  initial 
strain  method.  It  may  be  observed  that  it  differs  from  tangent  stiff- 
ness formulation  only  in  that  the  plastic  contributions  have  been  moved 

to  the  right-hand  side,  i.e.,  AF  - K Au,  leaving  K as  a constant 

~l>  «p  ~ «e 

global  3tiffness  matrix  that  requires  only  one  triangular ization.  To 
achieve  the  same  accuracy  as  the  tangent  stiffness  method,  it  is  neces- 
sary to  iterate  within  each  load  step  Many  of  the  calculations  are  the 
same  for  both  methods,  and  only  the  differences  are  emphasized  here. 

As  before,  It  is  assumed  all  quantities  are  known  at  time  step  n, 
and  the  objective  in  to  determine  the  quantities  at  time  step  n + 1. 
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1.  Determine  load  increment  AP;  let  u = u 

- ~n  -n-1 


n 

{-  a 

Estimate 

AF 

s y <bt  d 

B dv) 

Au 

~P 

at  m D 

V r 

% 

n ~n 

3. 

Solve  K 

Au 

AP  + AF  . If 

A A 

Au 

Au  , , go  to  Step  4;  other 

»e 

~n 

-P 

~n 

wise,  go 

to  Step 

2. 

4 

through  8. 

Same  as  tangent  stiffness 

method  * 

The  heart  of  the  initial  strain  method  is  embedded  in  the  iteration 

loop  within  the  load  step,  i.e.,  steps  2 and  3.  As  illustrated  in  the 

above  algorithm,  AF  is  reconstructed  on  each  iteration  based  on  the  new 
A ~P 

estimate  for  Au  . However,  there  is  no  reason  D cannot  also  be  changed 
~n  p 

simultaneously  on  each  iteration  to  reflect  abrupt  changes  in  transition 

zones  or  even  subtle  changes  in  D over  the  time  interval.  Herein  lies 

%P 

the  advantages  of  the  initial  strain  method.  To  wit,  after  the  initial 
triangular izat ion  of  Kg,  each  solution  merely  requires  modifying  the 
right-hand  side  and  performing  a back  substitution.  Moreover,  each 
estimate  of  AF^  can  simultaneously  consider  variations  in  D as  well  as 
Au  In  obtaining  a convergent  solution. 

This  concludes  the  development  and  implementation  of  plasticity 
theory.  In  Appendix  A some  common  loading  functions  are  presented 
together  with  some  fine  points  on  flow  laws  and  hardening  rules. 
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Chapter  3 


LINEAR  VISCOELASTICITY 

In  this  section  a general  tht ee-dimensional  isotropic  viscoelastic 
constitutive  formulation  will  be  developed  that  is  compatible  with  the 
finite  element  formulation  previously  presented.  To  this  end,  a brief 
review  of  viscoelastic  constitutive  theory  will  be  given,  beginning  with 
simple  one-dimensional  concepts  followed  by  a generalization  of  multi- 
dimensional stress-strain  states.  For  a comprehensive  introduction  into 
viscoelastic  theory,  the  reader  is  referred  to  Flugge  [19]  and  Christensen 
(20). 


BASIC  CONCEPTS 


Viscoelastic  materials  are  often  called  "memory"  materials,  that 
is,  the  stress  in  the  material  is  determined  not  only  by  the  current 
state  of  deformation,  but  also  by  all  past  deformation  states.  More- 
over, the  "memory"  exhibits  a fading  phenomenon  in  that  past  deformation 
states  influence  the  current  stress  state  to  a lesser  degree  than  do 
more  recent  deformation  states.  As  a consequence  of  this  memory  phenom- 
enon, viscoelastic  materials  dissipate  energy  during  deformat  ion;  thus, 
the  external  work  put  Into  the  system  cannot  be  completely  recovered. 
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ONE-DIMENSIONAL  MODELS 


The  above  characteristics  of  linear  viscoelasticity  can  be  exhibited 
by  one-dimensional  mechanical  models  composed  of  an  assemblage  of  springs 
and  dashpots,  where  the  springs  denote  a linear  relationship  between 
stress  and  deformation  while  dashpots  denote  a linear  relationship 
between  stress  and  rate  of  deformation.  Mechanical  models  provide  an 
insight  into  the  nature  of  viscoelastic  behavior  and  also  provide  a 
basis  for  defining  a differential  equation  relating  stress  and  strain. 

Differential  Equations 

Consider  the  one-dimensional  model  shown  below,  which  is  known  as 
the  standard  linear  solid. 


a = total  stress,  e = total  strain,  f.  = elastic  strain,  e = viscous 

e v 

strain,  = spring  constants,  and  u = dashpot  constant.  Using 

simple  concepts  of  equilibrium,  it  is  easy  to  deduce  ~ = in  the 


elastic  spring  and  c = 


in  the  parallel  assembly,  where  ( ) 


denotes  a time  derivative.  As  implied  by  the  model,  the  total  strain  is 
given  by  e = + e . Combining  these  three  relationships,  a differen- 

tial equation  relating  total  stress  to  total  strain  is  obtained  as: 
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0 + 


Ej  £ + PC 


(3-1) 


JL 

E1 


To  use  the  mechanical  model  as  a visual  aid,  it  is  noted  that  the 
dashpot  cannot  move  instantaneously.  Thus,  at  the  first  instance  of 
deformation  (from  t = 0 to  t = 0+)  only  the  spring  E^  can  deform; 
i.-»r.oef  o(0)  = e (0) . Furthermore,  if  the  load  is  maintained  for  an 
extended  period  of  time,  the  dashpot  will  finally  come  to  rest,  shifting 
all  the  load  to  spring  in  series  with  E^.  Hence,  at  t = «, 


o(») 


E1  E2 
E1  + E2 


e(°°) 


These  concepts  are  often  useful  in  interpreting  viscoelastic  responses. 

Equation  3-1  represents  a one-dimensional  viscoelastic  relationship 
for  a particular  model.  If  the  strain  were  prescribed  over  some  time 
interval  0 t <_  t , where  t is  the  current  time,  then  it  would  be  possi- 
ble to  solve  for  the  stress  response.  For  example,  suppose  the  strain 
is  input  by  a Heavyside  step  function,  such  that  c(t)  = h(t),  where 

c^  is  a constant  strain  magnitude  and  h(t)  is  the  unit  Heavyside  function 
(i.e. , h(t)  = 0 for  t < 0,  and  h(t)  = 1 for  t _>  0) , then  the  solution  of 
Equation  3-1  is: 


o(t) 


E1  E2 
E,  + E.. 


— exp[-t (E  + E9)/p] 
b2  1 


(3-2) 
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; where  the  initial  condition  o(0)  - E,  t was  used  to  determine  the 

1 o 

i constant  of  integration. 

! In  a manner  similar  to  the  above  example,  more  complicated  one- 

dimensional  models  can  be  constructed  by  adding  springs  and  dashpots  in 
series  and/or  parallel  that  will  always  result  in  a linear  differential 
equation  of  the  form: 

i 

i 

; p0  °0  + pi a + p2 H + •••  = % e + ‘‘i £ + ^2  *e'  + •••  (3“3) 

j 

I 

* where  p , p^  , p^,  ....  q^,  q^,  q2»  ...  are  constants  composed  of  the 

j parameters  of  the  springs  and  dashpots. 

Equation  3-3  may  be  considered  as  a general  linear  viscoelastic 

\ 

constitutive  relationship  for  a one-dimensional  problem.  However,  it  is 
not  a particularly  useful  form  for  solving  boundary  value  problems 
«,  because  both  stress  and  strain  are  expressed  in  terms  of  their  deriva- 

• ■ t j ves . 

Integral  Equations 

A more  useful  form  can  be  obtained  by  expressing  Equation  3-3  in  an 

' integral  form  in  the  following  manner.  First,  assume  the  strain  is 

prescribed  by  a Heavyside  step  function;  i.e. , t(t)  = h(t),  where 

is  some  constant  strain  amplitude  and  h(t)  = 0 for  t < 0 or  h(t)  = 1 for 

, t > 0.  With  this  prescription  for  strain,  the  right-hand  side  of  Fqua- 

. i t ion  3-3  reduces  to  the  constant  q e for  all  t > 0.  Therefore, 

o o 

Equation  3-3  represents  an  ordinary  linear  differential  equation  in 
terms  of  stress  fur  which  a solution  is  always  possible.  Symbolically, 
this  solution  can  be  expressed  as: 


j Ediki-  jilt, ,a ^Thiili £ jjj| 


Here  Y(t)  is  termed  the  relaxation  function  and  is  an  intrinsic  charac- 
teristic of  the  material  model  because  it  describes  the  nature  of  the 
stress  response  due  to  a prescribed  strain  equal  to  unity.  Implicit  in 
the  definition  of  a relaxation  function  is  Y(t)  = 0 when  t < 0;  i.e., 
there  is  no  stress  response  prior  to  application  of  strain.  As  an 
example,  the  relaxation  function  for  the  standard  linear  solid  is 
enclosed  in  braces  in  Equation  3-2. 

Next,  to  extend  the  solution  of  Equation  3-4  to  apply  to  an  arbi- 
trary strain  input,  linearity  and  superposition  are  exploited  as  follows. 
Consider  an  arbitrary  strain  input  e(t),  shown  in  Figure  12,  along  with 
a step-vise  approximation  to  this  curve  given  by 


l ! 


n 

e(t)  = ch(t)  + V Ac  h ( t - k At) 

k*l 

where  Ae^  is  an  increment  of  strain  applied  at  the  time  ",  = k At. 

Because  of  linearity  the  solution  given  by  Equation  3-4  applies  to 
each  strain  increment  corresponding  to  the  time  it  makes  its  contribu- 
tion; i.e.,  Ao^  = Y(t  - k A t ) A r ^ • Thus,  the  solution  for  stress  is 
given  by  the  superposition  of  solutions  as: 


o(t)  « Y(t)  e + Y (t  - x)  3--(^-d—  (3-6a) 

CJ  J 0 T 


For  simplicity  it  is  convenient  to  adopt  the  shorthand  notation  of 
convolution  algebra,  wherein  Equation  3-6a  is  written  equivalently  as: 

o(t)  = Y*  dc  (3-6b) 

where  * is  called  a convolution  operator. 

Equation  3-6a  or  3-6b  is  known  as  a heredity  integral  and  is 
physically  and  mathematically  equivalent  to  the  differential  form  given 
by  Equation  3-3.  However,  the  heredity  integral  provides  an  expression 
for  stress  without  stress  derivatives,  thereby  facilitating  displacement 
formulations  for  boundary  value  problems.  Moreover,  many  investigators 
assert  that  the  heredity  integral  is  the  proper  definition  for  the 
linear  viscoelastic  constitutive  law,  and,  thus,  there  is  no  need  to 
concoct  mechanical  models  or  their  associated  differential  equations. 

When  this  viewpoint  is  adopted,  the  relaxation  function  Y(t)  need  not  be 
considered  related  to  any  particular  mechanical  model,  but  rather  may  be 
considered  as  some  monotonically  decreasing  function  that  predicts  the 
stress  response  (relaxation)  which  occurs  from  an  imposed  unit  of  strain. 

In  the  remainder  of  the  study,  the  relaxation  function  will  be 
assumed  to  be  given  by  an  exponential  series,  i.e., 


where  E^,  E^,  E^,  ...  are  relaxation  moduli,  and  > ^ > 7^,  •••  are 
relaxation  tiroes.  Specification  of  these  non-negative  material  con- 
stants is  all  that  is  required  to  completely  define  the  relaxation 
function  which  in  turn  can  be  used  in  the  heredity  integral  to  form  a 
one-dime  isional  viscoelastic  constitutive  law. 

It  can  be  shown  that  the  exponential  series  form  of  th«  relaxation 
function  will  always  correspond  to  some  mechanical  model;  thus,  it  is 
left  to  the  discretion  of  the  reader  as  to  whether  or  not  he  will  accept 
the  relaxation  function  at  face  value,  as  given  in  Equation  3-7,  or 
interpret  the  relaxation  function  in  terms  of  a mechanical  model. 

To  summarize  the  discussion  thus  far,  Equation  3-6a  or  3-6b  repre- 
sents a general  one-dimensional  viscoelastic  constitutive  relationship, 
wherein  Y(t)  is  a relaxation  function  characterized  by  an  exponential 
series,  Equation  3-7.  The  parameters  of  the  exponential  series  may  be 
determined  directly  from  experimental  data  or  interpreted  from  mechanical 
models.  In  the  next  section,  multidimensional  stress-strain  models  will 
be  discussed. 

MULTIDIMENSIONAL  MODELS 

General  Viscoelastic  Constitutive  Law 

Extending  the  constitutive  relaxation  from  one  dimension  to  multi- 
dimensional stress  states  follows  reasoning  directly  analogous  o the 
generalized  Hook's  law  for  elastic  materials.  Namely,  each  component  of 
stress  is  coupled  in  some  fashion  to  various  components  of  the  strain 
vector  through  a constitutive  matrix,  i.e.,: 


o = D*  de 

m 


(3-8) 
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where  the  convolution  symbol,  *,  denotes  the  integral  relationships 
given  by  Equation  3-6b.  In  general,  the  matrix  D may  contain  21  inde- 
pendent relaxation  functions  for  describing  anisotropic  materials; 
however,  for  isotropic  materials,  only  two  independent  relaxation  func- 
tions are  required. 

In  this  writing,  only  the  isotropic  form  will  be  pursued.  It  is 
convenient  in  solving  boundary  value  problems  to  choose  the  isotropic 
relaxation  functions  as  the  response  due  to  bulk  and  shear  deformations. 
That  is,  the  bulk  relaxation  function,  K(t),  is  defined  as  the  hydrostatic 
stress  response  due  to  a prescribed  unit  of  volume  change  by  a Heavyside 
step  function.  Similarly,  the  shear  relaxation  function,  G(t),  is 
defined  as  the  shear  stress  response  due  to  a prescribed  unit  of  shear 
strain  by  a Heavyside  step  function. 

Thus,  for  isotropic  viscoelastic  materials,  Equation  3-8  has  the 
following  expanded  form,  where  K and  G represent  the  independent  relaxa- 
tion functions: 

K+(4G/3)  K-(2G/3)  K-(2C/3>  000 

K-(2G/3)  K+(4G/3)  K-<20/3)  000 

K-(2G/3)  K-(2G/3)  K+(4G/3)  0 0 0 

0 0 0 2G  0 0 

0 0 0 0 2G  0 

0 0 0 0 0 2G 


Equation  3-9  is  the  general  viscoelastic  constitutive  law  for 
isotropic  materials  and  has  a strong  resemblance  to  the  analogous  elastic 
constitutive  law.  However,  it  must  be  kept  in  mind  that  K and  G are 
functions  of  time,  and  the  convolution  operator,  *,  denotes  an  integral 
relationship.  For  example, 


By  integration  of  parts,  it  can  be  demonstrated  that  the  con- 
volution operator  is  communicative,  i.e.,  Y*  de  = e*  dY . Also,  the 
convolution  operator  is  linear  in  the  sense  that  if  Y = Y,r  + Y , then 
e*  dY  = e*  l + e*  dY  . Therefore,  Equation  3-9  may  be  uncoupled  in 
shear  and  bulk  and  written  as: 


o = tL(t*  dK)  + D„  (e*  dG) 


(3-10) 


This  may  be  written  in  expanded  form  as: 


Se  £(t>  - £k  /'  i(t> 

- Sa  f £<’> 


(3-11) 


where  DtJ 


1110  0 0 

1110  0 0 

1110  0 0 

0 0 0 0 0 .0 

0 0 0 0 0 0 

0 0 0 0 0 0 


(2/3) 


2-1-1000 
-1  2-1000 
-1-1  2000 
0 0 0 3 0 0 

0 0 0 0 3 0 

0 0 0 0 0 3 


K(°)  Dr 


G (o)  D_ 
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In  the  above,  D„  and  D„  are  constant  dimensionless  matrices  for 
K. 

bulk  and  shear,  respectively,  whereas  Dg  is  the  familiar  constant  elas- 
tic constitutive  matrix  representing  the  instantaneous  response  of  the 
viscoelastic  model. 

Viscoelastic  Law  With  Exponential  Series 


As  previously  discussed,  the  relaxation  functions  can  be  repre- 
sented by  an  exponential  series  with  no  serious  degradation  of  generality. 
Accordingly,  the  bulk  and  shear  relaxation  functions  are  taken  as: 


K(t) 


£ K exp ( — t / 6 ) 
i=l 


(3-12) 


G (t)  = G 


!£  G exp(-t/y  ) 
1=1 


(3-13) 


where  K^,  K , K_,  ...  Ku  and  G^,  G. , G»,  ...  Gv  are  relaxation  moduli 

1 K XL.  (.■» 


for  bulk  and  shear,  respectively,  and  6.,  ...  Bv  and  ,,, 

i Z lV  1 z 


Vjj  are  relaxation  times  for  bulk  and  shear,  respectively.  Determina- 
G 


tion  of  these  parameters  is  discussed  in  Appendix  B of  this  report  and 
elsewhere  [2]].  For  now  it  is  assumed  these  parameters  are  known,  and 
the  relaxation  functions  e completely  defined. 

Inserting  the  relaxation  functions,  i.e..  Equations  3-12  and  3-13, 
into  the  general  constitutive  law.  Equation  '3-11,  the  following  rela- 
tionship may  be  written: 
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£(t>  - £e  1U)  + £ Ki  hi  + £c  £-  Gi  Si 


(3-14) 


where 


Si 


e(0  -T-  exp[-(t 
Bi 


e(T) 


— exp [ -( t 

Yi 


T ) / B ± ] dT 


t ) / y t ] d i 


(3-15) 


(3-16) 


The  vector  sets,  k.  and  g.,  are  often  called  "hidden  coordinates" 

~i  “i 

or  "internal  variables"  (22,23)  for  bulk  and  shear,  respectively.  The 
motivation  for  introducing  these  internal  variables  will  become  evident 
in  subsequent  derivations.  It  will  be  shown  that  incremental  recursion 
relationships  for  k.  and  g can  be  used  to  great  advantage  in  circum- 
venting  the  need  for  storing  the  complete  history  of  deformation  [21,24] 

Incremental  Viscoelastic  Law 

In  anticipation  of  using  an  incremental  solution  procedure, 

the  constitutive  law  can  be  cast  in  incremental  form  by  defining  the 

current  time  increment  as  At  = t . - t and  Ac,  A •'  , Ak . , and  Ag,  as 

n+ 1 n - ~ ~i  Si 

A'  = c(t  ,.)  - "•  ( t ),  etc.  With  these  definitions.  Equation  3-14  can  be 
~ -»  n+1  ~ n ’ s 

directly  written  in  incremental  form  as: 


= 1) 


Ac 


Pv 


y k ak. 


Nc 


ug. 


(3-17) 


(iswffsrr 


where  1 - expC-At/y^) 

It  is  significant  to  note  that  Ak  and  Ag^  can  be  determined  by  jvist 
integrating  over  the  current  time  step,  t to  t j,  rather  than  the 
entire  time  history,  The  time  history  effect  is  given  by  the  recursive 
terms  rj  k^)  and  ^ g^^). 
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As  yet,  no  numerical  approximations  have  been  made.  However,  in 

order  to  evaluate  the  integral  in  Equations  3-18  and  3-19,  it  will  be 

assumed  the  strain  vector  varies  linearly  over  the  time  step.  That  is, 

in  the  time  interval  t < x < t + At , let: 

n — — n 


T - t 

c(t)  = e(t  ) + — ■ -"-Ae 

~ — n At  ~ 


Then  Equation  3-18  can  be  evaluated  as: 


Ak  = - r Ac  - r . [e(t  ) - k (t  ) 

~i  i ~ i ~ n -in 


(3-20) 


where 


ti  - 1 - exp (-At/ 8^) 


O 

, ' i - 

r . = 1 - tt  r 

i at  i 


In  an  identical  fashion,  an  incremental  recursion  relationship  for 
Ag  can  be  deduced  as: 


- qx 


- gi(tn)] 


(3-21) 


where 


= 1 - 


exp (-At / y ^ ) 


= 1 - 


'i  _ 
At  qi 
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Inserting  Equations  3-20  and  3-21  into  Equation  3-17  and  separating  the 
current  strain  increment  from  history  terms,  the  following  incremental 
constitutive  law  is  achieved: 


Aa  = (D  - D )Ae  - a 

— ' mV 


(3-22) 


where  D 


aK  £k  + UG  £g 


E k r 

i=l 


S ,, 

i=l  1 


o 

~v 


stress-history  influence  vector 


Note  that  D is  composed  of  the  constant,  dimensionless  matrices  D„  and 

nr  V m K 


0 , and  the  scalar  multiples  u and  a that  are  only  dependent  on  time- 


step  size  At  and  not  the  time  t.  Therefore,  D only  changes  its  value 

mt  V 


when  the  size  of  the  time  step  is  changed 

The  term  o will 

~v 

vector  and  is  given  by 


The  term  o will  be  called  the  viscoelastic  stress-history  influence 
~v  1 


!v(Ln)  = J.V  E Ki  - ~l(  tn) 


+ %G  £ Gi  i 1 tn)  ’ £i(tn} 


(3-23) 
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The  stress-history  influence  vector  accounts  for  the  influence  of 
all  past  deformation  states  on  the  current  stress  increment  given  in 


Equation  3-22.  It  should  be  observed  that  is  not  dependent  on  the 


current  time  interval  t to  t , , ; therefore,  a is  a known  quantity  at 

n n+1 ’ ~v  is 


the  beginning  of  each  time  step  and  can  be  treated  as  an  initial  stress. 


In  computational  practice  it  is  convenient  to  compute  at  the  end 


of  a time  step  in  preparation  for  the  next  step.  Summarized  below  are 

time  step  t to  t , : 

^ n n+1 


the  necessary  re 

lat ionships 

to  update  o after 

1. 

dhi  = 

- r 

Ac  - 

i — 

2. 

A2i  = 

- q 

, As 

- Si^n) 

3. 

*i<Cn+l 

) = 

W 

+ Ak. 

4. 

!i(Cn+l 

) = 

2i(tn} 

+ A2i 

5. 

£(tn+l> 

= 

E(tn) 

+ Ae 

6. 

£v(tn+l 

) = 

nk 

°K  £ 

i=i 

Ki  ‘ i 

)] 


+ £ Gi  WW  - ?i(t 


n+1 


n 


In  summary.  Equation  3-22  is  the  general  viscoelastic  constitutive 
relationship  to  be  incorporated  into  the  boundary  value  problem,  and  the 


above  relationships  provide  the  algorithm  for  updating  5 . In  the  next 


sections  the  constitutive  relationship  is  introduced  into  a finite 
element  formulation  and  a step-by-step  procedure  is  outlined  for  the 
solut ion. 


•i 

'im 

! 


3. 

i 
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VISCOELASTIC  FINITE  ELEMENT  FORMULATION 


For  convenience,  the  general  equilibrium  equations  previously 
derived  in  the  finite  element  development  are  repeated  here; 


Ao  dv  = AP 


(3-24) 


where  it  will  be  recalled  that  B is  the  strain-to-nodal  point  displace- 
ment matrix,  AP  is  the  external  load  increment,  and  the  symbol  & 
denotes  the  ordered  summation  of  each  finite  element  volume  integration 
into  its  correct  location  of  the  global  stiffness  matrix.  Replacing  Ao 
by  the  viscoelastic  constitutive  law.  Equation  3-22,  and  using  the 
strain-to-node  displacement  relationship,  A e = B Au,  Equation  3-24  may 
be  written  as: 


(K  - K )Au  = AP  + F 
»e  »v  ~ ~v 


(3-25) 


where 


K = / BT  D b dv 

~e  J m „e  . 

V 

K = C / B1  D B dv 

-v  L-!  j - -V 
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F 

~v 


(3-28) 


= 2 / £ 2V  dv 

v 


Equation  3-25  represents  the  familiar  set  of  linear  algebraic 

equations  common  to  step-by-step  methods.  is  the  constant,  elastic, 

global  stiffness  matrix,  and  is  the  global  viscous  stiffness  matrix 

which  is  independent  of  time,  t,  and  is  only  dependent  on  the  current 

size  of  the  time  step,  At.  F^  represents  the  viscoelastic  force  history 

vector  which  is  independent  of  the  current  time  interval,  t to  t ,, 

n n+1 

and,  therefore,  is  known  at  the  beginning  of  each  time  step. 

The  solution  procedure  for  Equation  3-25  can  be  handled  in  a 
straightforward  manner.  However,  caution  must  be  followed  with  regard 
to  the  calculation  of  F ^ . Note  the  volume  integral  in  Equation  3-28 
requires  that  the  spatial  distribution  of  be  known  within  each  ele- 
ment. Unfortunately,  it  is  not  known  in  a continuous  fashion.  There- 
fore, in  order  to  perform  the  volume  integration,  several  approaches  are 
possible.  First,  it  may  be  assumed  that  o^  is  constant  within  the 

element.  Thus,  o is  determined  at  the  element  center.  Second,  o mav 
~v  ~v 

be  computed  at  the  numerical  integration  points  (Gauss  points) , thereby 
facilitating  a direct  numerical  integration.  Although  this  second 
approach  is  accurate,  it  may  require  excessive  auxiliary  storage  because 
the  calculation  of  3 requires  the  storage  of  all  the  internal  variables 
k^  and  g for  each  integration  point  (see  Equation  3-23). 

A third  and  last  approach  is  the  most  consistent.  Basically,  the 
idea  is  to  develop  an  updating  procedure  for  F^  based  on  a recursion 
scheme.  To  this  end,  Equations  3-23  and  3-28  are  combined  to  express  F 
in  the  following  manner: 


F 

~v 


K 

L 

i=l 


+ 


(3-29) 
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where 


M 

t ! 


V = Ki  7i  If  - ii(tn) 


]dv 


(3-30) 


*Gj  = Gi  qi  /f  °G^(tn}  " h(tn 


) 3d v 


(3-31) 


Tn  Che  above  equations,  h^  and  represent  viscoelastic  force  history 
vectors  for  bulk  and  shear  within  each  element. 


As  expressed  in  Equations  3-29,  3-30,  and  3-31,  the  calculation  of 
and  hg^  and  then  F^  does  not  provide  any  advantage  over  the  calcula- 
tion for  F as  given  previously  in  Equation  3-28.  However,  if  Equations 


3-30  and  3-31  are  written  in  Incremental  form  (that  is,  Ain-  = hi<  (t  ) - 

n 


hv  (t  ,)  and  Ahr,  = hr,.(t  ) - hc.Ct  ,),  if  Ac  is  replaced  by  B Au,  and 
n-1  n -i  n-1  ~ - 


if  Ak  and  Ag,  are  replaced  by  their  recursion  relationships  given  by 

• i a«i 


Equations  3-20  and  3-21,  then  the  following  recursive  relationships  are 
found  for  Ahv-  and  Ah,~  . 

<w  N J 


\ “ Ki  V1  + ri>SK  ^ 


~K . 


(3-32) 


AtGi  = Gi  qi(1’+  VSc 


qi  ~G 


(3-33) 


ess 


f T 

where  Sv  = / B D B dv , dimensionless  bulk  stiffn 

Mv  K.  v Mr  Mi  K % 

V 

f T 

S_  = ] B D B dv , dimensionless  shear  stiffness 


t 

[ 


At  the  end  of  each  time  step.  Equations  3-32  and  3-33  provide  a 
simple  and  consistent  method  of  determining  the  force  history  increments 
Ahe  and  A he,  for  each  element  by  making  use  of  the  known  displacement 
increment  Au  and  the  accumulated  force  history  vectors,  h^  and  h^ . . 

It  should  be  observed  that  the  volume  integrations  are  consistently 

defined  bv  the  dimensionless  element  stiffness  matrices,  S,,  and  S_. 

•>  K ‘•G 

furthermore,  S..  and  S..  are  independent  of  time  and  time  step;  therefore, 

mG 

they  need  only  be  computed  once  for  each  element  and  stored. 

The  foregoing  viscoelastic  finite  element  formulation  is  summarized 

in  the  following  step-by-step  solution  procedure.  Equation  3-25  is  the 

governing  equilibrium  equation  to  be  solved  at  each  time  step.  To 

present  the  solution  strategy,  it  is  assumed  that  ail  quantities  have 

been  calculated  at  time  t , i.e.,  u(t  ),  h,,  (t  ),  and  hr  (t  ).  The 

n ~ n ~|'|  n n 

obiective  is  to  find  the  increments  Au,  Alw  , and  Aho  . from  time  t to 

~ n 

tn+l ' 

1.  Form  new  load  increment  AP. 

2.  Assemble  force  history  vector,  , from  Equation  3-29. 

3.  Assemble  total  stiffness  matrix,  K - K , from  Equations  3-26 

<xe  *>v  ^ 

and  3-27  (only  required  if  time-step  size  changes). 

4.  Solve  (K  - Y.  )Au  = AP  + F for  Au. 

_e  «v  ~ ~ -v  ~ 

5.  Evaluate  Ah  and  Ah  from  Equations  3-32  and  3-33. 

. — K , 

1 1 


6.  Update  all  quantities:  u(tn+1),  h(,  (tn+^),  h (tn  ). 

~ ’i  i 

7.  Print  results,  and  return  to  Step  1.  (Note  stresses  may 
be  .-dlculated  from  Equations  3-22  and  3-23.) 

To  start  the  above  algorithm,  it  is  assumed  the  body  is  undeformed 

0.  For  the 


prior  to  loading  so  that  initially  u(o)  = (o)  = h^.(o) 

initial  instantaneous  elastic  solution  over  the  time  interval  t - 0 


to 


0 , we  have  At  = 0,  which  implies  K = 0,  Ah  = Ah  = 0,  and  F = 0. 

m V K. . . *»V 

1 1 


L 
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Therefore,  the  first  solution  at  t = 0+  is  obtained  by  solving  the 
elastic  system,  Ke  Au  = AP , where  AP  is  the  instantaneous  load  applied 
at  t = 0.  With  this  starting  procedure,  the  above  algorithm  can  be  used 
for  each  succeeding  step. 

As  previously  mentioned,  the  viscoelastic  stiffness  matrix, 
only  changes  its  value  when  the  size  of  the  time  step,  At,  changes.  If 
the  time  step  is  kept  constant,  computation  time  (steps  3 and  4)  may  be 
considerably  reduced  by  reusing  the  total  tr iangularized  stiffness 
matrix  to  modify  each  new  right-hand-side  vector,  AP  + F , and  perform- 
ing a simple  back  substitution  to  determine  Au. 
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Chapcer  4 
VISCOPLASTICITY 

i 

GENERAL 

The  general  theory  of  viscoplasticity  is  discussed  in  Appendix  C. 

In  this  section,  a particular  viscoplastic  model  called  "viscoelastic- 
: plastic"  is  developed  in  detail.  The  model  is  a member  of  a family  of 

combo-viscoplastic  models  introduced  in  Appendix  C. 

The  viscoelastic-plastic  model  is  characterized  by  a linear 
viscoelastic  model  within  the  yield  surface  and  the  combined  visco- 
elastic and  plastic  response  on  the  yield  surface.  Figure  13  portrays 
a one-dimensional  representation  of  the  viscoelastic-plastic  model.  As 
suggested  by  the  one-dimensional  model,  plastic  deformation  is  not 
retarded  by  viscous  components;  consequently,  plastic  deformation  may 
occur  instantaneously  as  in  classical  plasticity.  Accoruingly,  the 
yield  function  is  restricted  to  the  same  rules  of  classical  plasticity. 
The  viscoelastic-plastic  model  is  developed  in  detail  in  the  next  section 
for  a general  multidimensional  stress-strain  state  wherein  previously 
derived  relationships  for  plasticity  and  viscoelasticity  are  employed. 


VISCOELASTIC-PLASTIC  CONSTITUTIVE  DEVELOPMENT 

In  this  section,  a general  constitutive  relationship  is  developed 
for  the  class  of  "viscoelastic-plastic"  materials  described  above.  The 
fundamental  assumption  for  viscoelastic-plastic  materials  is  that  the 
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total  stress  is  persistently  related  to  the  viscoelastic  portion  of  the 
total  strain  regardless  of  the  amount  of  plastic  straining,  i.e., 


a = D*  de 
~ « ~ve 


(4-1) 


where  D*  devg  denotes  the  convolution  constitutive  relationship  pre- 
viously developed  in  Equations  3-8  through  3-11. 

In  accordance  with  the  above  assumption,  it  is  further  assumed  that 
the  total  strain  can  be  decomposed  into  viscoelastic  and  plastic  compo- 
nents , i.e., 


e = 


e 

~ve 
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(4-2) 


where  e is  the  viscoelastic  strain  vector,  and  e is  the  plastic 
~ve  «*p 

strain  vector.  Note  the  assumptions  embodied  in  Equations  4-1  and  4-2 
are  directly  analogous  to  the  assumptions  made  in  the  elastic-plastic 
formulation,  wherein  it  is  assumed  a = D e and  e = e + e . 

As  always,  the  objective  of  this  constitutive  development  is  to 
determine  a relationship  between  the  current  stress  increment,  Aa , and 
the  strain  increment,  Ac.  To  this  end,  the  development  begins  with  the 
incremental  equivalent  of  Equation  4-1  given  by  Equation  3-22  and  is 
repeated  here  for  convenience: 
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(4-5) 


In  Equation  4-3,  the  incremental  quantities  imply  the  increments  from 

time  t to  t , , . The  stress-history  influence  vector,  a , is  evaluated 
n n+1  J ’ ~v 

at  time  t and,  therefore,  can  be  considered  as  a known  initial  stress, 
n 

All  of  the  quantities  on  the  right-hand  side  of  Equations  4-4  and  4-5 
were  previously  defined  and  discussed  in  the  viscoelastic  formulation. 
The  significant  point  to  be  borne  in  mind  is  that  all  references  to 
strains  in  Equations  4-3  and  4-5  are  with  respect  to  the  viscoelastic 
strain  and  not  the  total  strain. 

To  achieve  the  desired  constitutive  form,  A e may  be  replaced  by 

Ac  - Ac  in  accordance  with  Equation  4-3;  thus.  Equation  4-3  becomes: 

~ ~P 


Acs  = D (Ac  - Ac  ) — o 

- «ve  ~ ~p  «*v 


(4-6) 
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Except  for  the  terra  5^,  Equation  4-6  Is  similar  to  Equation  2-30  in  the 
elastic-plastic  development.  Therefore,  at  this  point,  the  development 
parallels  the  plastic  formulation,  and  it  is  assumed  the  plastic  flow 
law  (Equation  2-26)  remains  valid,  i.e.. 


B A£S  a 

H*  B 


where  n is  the  outward  normal  of  the  yield  surface,  m is  the  direction 
of  plastic  straining  (m  = n implies  associative  law;  m 4 n implies 

•»  mm  am 

nonassociative  law),  and  H'  is  the  hardening  coefficient.  Furthermore, 
it  is  assumed  all  the  previously  established  rules  for  plastic  deforma- 
tion still  hold.  That  is,  yielding  can  only  occur  for  states  of  stress 
on  the  yield  surface  with  some  component  of  the  stress  increment  colinear 
with  n,  and  the  methodology  of  tracking  the  universal  yield  surface 
remains  unaltered. 

Inserting  the  flow  rule,  Equation  4-7,  into  Equation  4-6  and  taking 

aT 

the  inner  product  with  respect  to  the  outward  normal,  n , the  scalar 

aT  ~ 

quantity  n A a is  determined  as: 
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(4-8) 


Inserting  Equation  4-8  back  into  Equation  4-7,  the  plastic  strain  incre- 
ment is  related  to  the  total  strain  increment  by: 


At  - C (D  At  - o ) 
-p  «pv  *ve  - ~v 


(4-9) 
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Lastly,  returning  to  Equation  4-6  and  replacing  Ac  by  Equation  4-9,  the 
desired  viscoelastic-plastic  constitutive  relationship  is  achieved: 


Ao  = (D  - D )Ae  - 5 

~ »ve  *vp  ~ ~vp 


(4-10) 


where 


D = D T C D 
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o = (D  - D )D  1 a 
~vp  *ve  «%vp  «ve  ~v 


(4-12) 


In  the  above,  D is  a viscoplastic  material  matrix,  and  o is  the 
*>vp  ~vp 

v i scop  Las tic  stress-historv  influence  vector  and  is  related  to  o 

~ v 

through  Equation  4-12.  Note  o reduces  to  c whenever  D =0. 

~vp  ~v  »vp  » 

In  utilizing  Equation  4-10,  o is  treated  as  a known  initial 

~vp 

stress  during  each  time  interval.  At  the  end  of  each  time  step  o is 

~vp 

updated  in  preparation  for  the  next  step.  Summarized  below  are  the 
necessary  relationships  to  update  a after  A r and  Aa  have  been  deter- 
mined in  the  time  interval  t to  t 
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In  summary,  Equation  4-10  is  t he  general  viscoelastic-plastic 
incremental  constitutive  r»'l;  tionsbip  useable  in  any  boundary  value 
i ■•’■mulation.  The  above  sequence  of  steps  provide  the  algorithm  for 
updating  5 at  the  end  of  eacli  time  step.  Lastly,  it  can  be  observed 
th't  Equation  4-10  reduces  to  the  linear  viscoelastic  model  (Equation 

4-3)  whenever  D =0  (i.e.,  C = 0).  Alternatively,  if  viscous  compu 

~vp  « „p 

nents  arc  zero  (i.e.,  D =1)  and  n =0),  then  Equation  4-10  reduces 

•ve  me  ~v 

to  the  elastic-plastic  model. 

In  the  next  section,  Equation  4-0  is  incorporated  into  a finite 
element  formulation,  and  a step-by-step  solution  procedure  is  pi  •■.■nte 
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FINITE  ELEMENT  VISCOELASTIC-PLASTIC  FORMULATION 

beginning  with  the  general  finite  element  equilibrium  equations 
(from  Equation  1-6): 


Ao  dv  = AP 


(4-13) 


where  B is  the  strain-to-nodal  point  displacement  matrix  (i.e., 

Ac  = b*Au),  AP  is  Che  external  load  increment,  and  § implies  the  ordered 
summation  of  the  elements. 

Replacing  Ao  by  the  viscoelastic-plastic  constitutive  law,  Equation 
4-10,  and  using  Ac  * B Au,  Equation  4-13  may  be  written  as: 

* tm  KH 
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In  the  above,  Kve  Is  the  global  viscoelastic  stiffness  matrix  and 

could  be  written  as  K = K + K as  was  done  in  the  viscoelastic  formu- 

«ve  %v  «e 

lation.  Recall  that  K remains  constant  if  the  time  step  is  constant. 

K is  the  global  plastic  stiffness  (reduction)  matrix  and  is  similarly 

defined  in  the  plastic  formulation.  The  vector  F will  be  called  the 

~vp 

viscoelastic-plastic  force  history  vector  and  represents  the  interaction 

of  viscous  and  plastic  response.  If  D * 0,  F reduces  to  the  pre- 

■»  v p -»vp  ^ 

viously  defined  viscoelastic  force  history  vector,  Fv» 

The  governing  equations,  as  shown  in  Equation  4-14,  are  in  the 

proper  form  for  a direct  solution  by  standard  elimina'ion  methods. 

Alternatively,  the  plastic  portion,  K Au,  could  be  taken  to  the  right- 

hand  side  and  treated  as  an  unknown  force.  These  solution  strategies 

will  be  detailed  later;  for  now  the  attention  is  focused  on  computing 

the  force  history  vector,  F 

3 -vp 

The  calculation  of  Fv  presents  the  same  problem  as  did  F^  in  the 

viscoelastic  development;  namely.  Equation  4-17  requires  the  volume 

Integration  of  o^;  however,  the  spatial  distribution  of  o is  not 

known  in  a continuous  fashion.  There  are  two  basic  options  available  to 

overcome  this  difficulty.  One  method  is  to  calculate  o at  a suffi- 

~vp 

clent  number  of  points  to  afford  a proper  representation  in  a numerical 
integration  scheme.  For  example.  In  Gaussian  quadrature,  the  calcula- 
tion of  at  the  Gauss  points  would  be  appropriate.  The  disadvantage 
of  this  method  is  chat  each  point  at  which  a is  calculated  requires 
the  storage  of  all  the  Internal  variables  k and  g,. 

The  alternative  method  is  to  use  element  force  history  vectors 

(e.g.,  hi;  and  hr  ) as  was  done  in  the  viscoelastic  formulation.  How- 
* i **  i 

ever,  in  this  case,  the  development  is  not  as  clean  due  to  the  interaction 

of  viscous  and  plastic  responses.  As  a result,  a portion  of  F is 

**  * P 

dependent  on  pLastic  strain  which,  in  turn,  must  be  evaluated  at  points 

within  the  element  and  numerically  integrated.  Consequently,  as  long  as 

it  is  necessary  to  calculate  plastic  strains  at  the  integration  points, 

it  is  only  a little  more  work  to  obtain  o at  Gauss  points  and  obtain 

„vp 

F by  the  first  method. 

-vp 
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SOLUTION  ALGORITHM 


I 

! 

i 


* 
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The  solution  procedures  developed  in  the  elastic-plastic  and 
viscoelastic  sections  contain  most  of  the  concepts  for  solving  the 
viscoelastic-plastic  problem.  Here,  only  the  highlights  of  the  initial 
strain  method  will  be  discussed. 

From  Equation  A-1A,  the  proper  form  for  the  initial  strain  method 
is  given  as: 


K Au  = AP  + AF  + F ( A-18) 
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Recall  that  remains  constant  as  long  as  the  size  of  the  time  step 

remains  constant:  thus,  K only  needs  to  be  triangular ized  when  the 

<ve 

time  step  size  changes.  To  establish  the  algorithm,  it  is  assumed  all 
quantities  are  known  at  time  t^,  and  the  objective  is  to  determine  the 
quantities  at  time  cn+^* 

1.  Form  new  load  increment  AP,  and  estimate  Au , = Au  , 

~ ~i  ~n-l 


2.  Calculate  AF 

-P 


E(  J RT  D B dv\  Au 
^ J m «vp  . J -i 


3.  Compu 


te  F = S / fiT  ° 

~Vp  i—J  •/  w ~Vp 


dv,  where  o is  given  bv  Equa- 
~vp  ^ 


tion  4—12 
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4.  Solve  K Au.  = AP  + AF  + F for  Au , 

-ve  ~i  ~ ~p  ~vp  ~i 

5.  If  Au^  = Au^  go  to  Step  6;  otherwise,  return  to  Step  2 
for  iteration. 

6.  Monitor  each  element  to  determine  if  it  is  in  plastic  region. 

If  so,  update  yield  surface  parameters  and  flow  law. 

7.  Print  desired  results,  and  return  to  Step  1. 

The  above  algorithm  implies  D and  o remain  constant  within  the  time 

mV p ~Vp 

step;  however,  these  quantities  could  be  modified  in  the  iteration  loop 
to  account  for  transitions  from  the  viscoelastic  range  to  the  plastic 
range  as  was  discussed  in  the  section  on  plasticity. 

This  concludes  the  viscoelastic-plastic  development.  It  is  felt 
this  model  incorporates  all  of  the  significant  aspects  of  viscoelasticity 
and  plasticity  into  a unique  model.  Furthermore,  the  model  directly 
reduces  to  plasticity  or  viscoelasticity  by  appropriate  designation  of 
the  parameters. 
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Chapter  5 


II ENTIFICATION  AND  APPLICATION  OF  VISCOELASTIC-PLASTIC  MODEL 


In  this  section,  the  viscoelastic-plastic  model  is  discussed  from  a 
conceptual  viewpoint  to  illustrate  the  types  of  material  behavior  that 
can  be  represented  and  the  methods  for  determining  parameters.  To 
supplement  the  discussion,  the  model  is  compared  with  experimental  data 
from  sea-ice  and  plexiglass  test  specimens. 

GENERAL  CHARACTERISTICS 


A conceptual  representation  of  the  viscoelastic-plastic  model  is 
shown  in  Figure  14  where  the  symbols  (?) , © , and  © represent  the 
plastic,  viscous,  and  elastic  parts  of  the  model.  As  suggested  by  the 
figure  the  total  strain,  e,  is  the  sum  of  the  parts  (e  = e + e + e ). 

* nf  * ^ ryj  g 

Also,  the  figure  illustrates  that  the  applied  stress,  o,  is  transmitted 
from  (?)  to  © to  (?)  instantaneously  at  full  value,  because  the  com- 
ponents are  all  in  series.  Indeed,  these  characteristics  were  used  in 
the  previous  section  for  developing  the  general  constitutive  model. 

As  a side  comment,  other  combinations  of  models  with  different 
characteristics  can  be  developed  by  various  arrangements  of  the  basic 
components.  This  is  discussed  in  Appendix  C. 

To  illustrate  the  viscoelastic-plastic  model  shown  in  Figure  14, 
consider  the  strain  response  due  to  a constant  stress,  o = (l/2)c^  , 
where  is  the  linear  limit.  Since  o is  within  the  yield  surface,  only 
elastic  and  viscous  deformation  will  occur.  This  is  illustrated  by  the 
bottom  curve  in  Figure  15.  Note  the  initial  response  is  the  instantaneous 
elastic  strain,  and  all  additional  strain  accumulation  is  viscous. 


Figure  14.  Conceptualization  of  viscoelastic-plastic  model. 
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Figure  15.  Typical  strain  histories  for  various 
levels  of  stress. 
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Next,  consider  the  same  test  with  a = (slightly  less).  Again, a is 
still  within  the  yield  surface.  Therefore,  the  strain  response  is 
double  the  previous  strain  response  as  shown  by  the  next  curve.  Of 
course,  this  result  is  due  to  the  superposition  principle  of  linear 
viscoelasticity. 

Now  consider  a = (3/2) oT . This  time  the  yield  stress  is  exceeded, 
triggering  a plastic  response.  Like  the  elastic  strain,  the  plastic 
strain  occurs  instantaneously  and,  thereafter,  remains  constant  because 
it  only  responds  to  stress  changes.  The  magnitude  of  the  plastic  strain 
is  dependent  on  the  hardening  function,  shape  of  the  yield  surface,  and 
level  of  stress.  Figure  15  shows  a representative  strain  history  for 
this  example.  As  the  stress  level  is  increased  on  each  subsequent  test, 
the  strain  magnitudes  increase  in  a nonlinear  fashion  until  the  hardening 
function  becomes  zero,  whereafter  unrestrained  plastic  deformation 
occurs.  Again,  Figure  15  illustrates  these  concepts. 

A significant  characteristic  of  this  model  is  that  for  states  of 
constant  stress,  only  the  viscous  strain  changes  with  time  so  that  the 
strain  rate  is  independent  of  plastic  deformation  and  is  linearly  related 
to  stress  level  like  an  ordinary  viscoelastic  model. 

The  above  insights  are  useful  for  determining  the  parameters  of  the 
model  discussed  next. 


MODEL  FITTING  TECHNIQUES 

The  beauty  of  the  viscoelastic-plastic  model  presented  herein  is 
that  it  is  a combination  of  two  well-known  constitutive  theories:,  visco- 
elasticity and  plasticity.  Accordingly,  the  material  information  gathered 
over  the  years  for  the  individual  theories  can  be  used  directly  for  the 
combined  model. 

The  first  step  in  establishing  a viscoelastic-plastic  model  repre- 
sentative of  a particular  material  is  to  determine  the  viscoelastic 
portion  of  the  model  (i.e.,  (e)  + (v)  of  Figure  14).  If  the  material  is 


i 
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well  studied,  it  is  likel-r  that  viscoelastic  parameters  can  be  found  in 
the  engineering  literature.  Ocherv;ise,  Appendix  B provides  a detailed 
guideline  for  converting  test  data  into  viscoelastic  parameters. 

In  this  case,  the  most  important  consideration  is  to  use  test  data 
with  state  variables  (i.e.,  load  rates,  temperature,  moisture,  etc.)  as 
close  as  possible  to  the  conditions  of  the  problem  being  studied. 

Since  the  viscoelastic  model  represents  the  linear  portion  of  the 
response,  it  should  be  based  on  test  data  at  low  stress  levels.  Or, 
more  ideally,  it  should  be  based  on  an  unloading  response  history  since 
unloading  is  assumed  to  be  linear. 

The  plastic  portion  of  the  model  may  be  simple  or  as  complicated  as 
desired.  In  general,  a yield  condition,  flow  rule,  surface  (hardening) 
rule,  and  hardening  function  must  be  determined.  For  many  materials 
these  are  presented  in  the  engineering  literature.  Also,  Appendix  A 
provides  a detailed  discussion  on  various  plasticity  models. 

If  very  little  is  known  about  the  material  being  considered,  the 
plasticity  model  should  be  made  as  simple  as  possible;  i.e.,  assume  a 
simple  yield  condition  (s-ich  as  Drucker-Prager  or  Von  Mises),  assume  an 
associative  flow  rule,  assume  an  isotropic  hardening  surface,  and  calcu- 
late a hardening  function  directly  from  available  test  data. 

To  illustrate  the  viscoelastic-plastic  model  fitting  procedure, 
consider  the  idealized  strain  history  data  curves  that  were  discussed  in 
Figure  15.  The  first  step  is  to  determine  the  viscoelastic  creep  func- 
tion J(t),  satisfying  the  relationship  ?(t)  = J(t)c,  (c  < c ) . In  other 
words,  by  inspecting  the  data,  one  recognizes  the  lower  two  curves  obey 
the  linear  viscoelastic  superposition  principle  so  that  J(t)  may  be 
determined  from  either  one  of  these  curves.  The  detailed  process  o; 
determining  J't)  is  given  in  Appendix  B,  where  J(t)  is  of  the  form: 


N 

= A + B t + T C , ( 1 


J(t) 


e 


(5-1) 


and  the  unknown  parameters  to  he  determined  are  A,  B,  C ^ , C2>  C^, 

Bl*  6?.’  ' ’ ‘ BN' 

To  identify  a complete  viscoelastic  model,  additional  material  data 
are  required  (e.g.,  shear  or  bulk  function)  as  discussed  in  Appendix  B. 
However,  if  no  other  data  are  available,  some  constant  value  of  Poisson's 
ratio  can  be  assumed.  Methods  of  inverting  creep  functions  to  relaxation 
functions  and  converting  uniaxial  functions  to  bulk  and  shear  functions 
are  also  presented  in  Appendix  B. 

After  the  viscoelastic  parameters  are  determined,  the  plastic 
portion  of  the  model  is  examined  by  constructing  a plot  of  stress  versus 
plastic  strain.  To  this  end,  consider  any  particular  time  t*  in  Figure 
15.  By  definition,  the  plastic  strain  is 


tp  = G(C*>  - £vgU*) 


where  c(t*)  is  the  strain  data  at  some  stress  level  o,  and  e (t*)  * 

ve 

Therefore,  for  any  time  t*  a plot  of  stress  versus  plastic 
strain  may  be  constructed  as  shown  in  Figure  16.  If  the  material  is 
ideally  viscoelastic-plastic,  this  plot  will  be  identical  for  any  choice 
of  t*.  However,  in  practice,  it  is  prudent  to  consider  several  different 
values  of  t*  to  get  an  average  overtime.  In  so  doing,  one  may  discover 
the  plastic  strain  is  highly  time-dependent  so  that  a viscoelastic- 
plastic  model  is  not  suitable.  In  this  case,  consult  Appendix  C for 
other  combo-viscoplastic  models. 

Assuming  the  plastic  strain  remains  reasonably  constant  for  each 
stress  level,  the  next  step  is  to  choose  a plastic  yield  condition. 
Generally,  the  Drucker-Frager  yield  condition  (see  Appendix  A)  is  suffi- 
cient to  characterize  material  behavior  from  one-dimensional  tests. 
Assuming  isotropic  hardening,  this  condition  is: 
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where  a is  the  average  hydrostatic  stress,  and  is  the  second  devia- 
m l 

toric  stress  invarier.t.  The  yield  surface  may  be  visualized  in  principal 
stress  space  as  a cone  whose  centerline  is  along  the  hydrostat  (rr  = 

°2  = °q)*  The  yield  parameter,  k,  is  the  cone  radius  measured  from  the 
origin,  and  the  parameter  -i  measures  the  rate  of  increase  of  the  cone 
radius  in  moving  along  the  compression  hydrostat.  If  a = 0,  the  cone 
becomes  a right  cylinder  and  the  Drucker-Prager  condition  reduces  to 
the  Von  Mises  yield  condition. 

For  the  problem  under  consideration,  the  data  in  Figure  16  contain 
two  compression  tests  that  exceed  the  elastic  limit.  However,  because 
these  Lests  have  identical  stress  paths  and  intercept  t he  same  point  on 
the  yield  surface,  they  do  not  provide  sufficient  information  to  deter- 
mine both  k and  n independently.  If  tension  yield  tests  are  conducted 
in  addition  to  the  compression  yield  tests  (or,  if  lateral  plastic 
strains  are  measured),  then  both  k and  a can  be  determined. 

Assuming  for  now  only  compression  yield  data  are  available  (which 
is  often  the  case),  then  a must  be  pre-selected  as  some  constant,  say 
a = (most  probably  u = 0).  From  Figure  16,  initial  plastic  yielding 
occurs  at  r.  = o.  so  that  from  Equation  5-?  the  initial  value  of  k is: 


k — — — o + 

o 3 L 


(5-3) 


where  < 0 is  the  initial  compressive  yield  stress.  During  plastit 
yielding,  k increases  such  that  it  continuously  satisfies  the  yield 
cond  Ltion. 


Plastic  Strain,  Cp 


Figure  16.  Typical  plastic  hardening  curve. 


Having  established  a yield  condition  and  its  parameters,  one  next 
assumes  an  associative  flow  rule  and  determines  the  hardening  function. 
From  Equation  2-24,  the  associative  flow  rule  for  the  plastic  strain 
increment  gives  the  hardening  function  as: 


H(W  ) 
P 


(5-4) 


where  is  the  component  of  the  yield  surface  unit-normal  collateral 
with  the  loading  stress,  and  is  given  by: 


n 
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(5-5) 


If  a =0  (i.e.,  Von  Mises  condition),  the  above  reduces  to  n = (2/  V6) 
o l 

Sgn  a. 

Using  Equations  5-4  with  n^  given  by  Equation  5-5,  and  Ao/ given 
by  the  slopes  in  Figure  16,  the  hardening  function  can  he  evaluated  as  a 
function  of  plastic  work 


W 

P 


dL 


p 


(i.e.,  area  under  a versus 

for  the  case  a =0. 

o 


c 

P 


curve)  . 


This  is  illustrated  in  Figure  17 
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The  preceding  illustrates  the  basic  strategy  of  identifying  the 
parameters  of  the  viscoelastic-plastic  model.  Appendixes  A and  B pro- 
vide further  insight  and  details.  In  the  next  section  some  examples  are 
given . 


COMPARISON  WITH  EXPERIMENTAL  DATA 

In  this  section,  the  versatility  (and  limitations)  of  the 
viscoelastic-plastic  model  is  demonstrated  by  comparing  it  to  experi- 
mental strain  data  from  compression  test  specimens  of  Plexiglas  and  sea- 

ice  . 

The  Plexiglas  creep  data  were  generated  by  Marin,  Pao,  and  Cuff 
[271,  and  are  shown  by  discrete  data  points  in  Figure  18  for  three 
compressive  stress  levels,  c = 2,500,  3,100,  and  4,000  psi.  Using  the 
model  fitting  procedure  in  t Vie  previous  section  (including  Appendixes  A 
and  B),  the  viscoelastic-plastic  model  expression  for  creep  strain  is: 


,(t)  = —-r  [2.57  + 0 . 000417 1 + 0.389(1  - e °'3t)]  + • 

io6  P 


with  l =0,  when 

P <• 


_ , when  - oL 
P P 


!t  is  assumed  = * . 1 on  (initial  linear  limit),  ar.u  EC  ' = 50°, 000 

I . P p 

psi  is  simplified  to  a constant.  The  time  unit  is  hours. 


The  above  expression  is  plotted  as  a solid  line  together  with 
the  experimental  data  points  for  all  three  stress  levels.  For  the  two 
lower  stress  levels,  the  model  is  linear  viscoelastic  (i.e. , no  plastic 
response)  and  compares  well  with  the  data  points.  For  the  highest 
stress  level  (4,000  psi)  a plastic  strain  is  generated  that  combines 
with  the  viscoelastic  strain.  This  gives  a better  fit  of  experimental 
strain  than  just  the  linear  viscoelastic  model  as  shown  bv  the  dashed 
line. 


The  above  example  illustrates  a predominantly  elastic  response,  a 
moderate  amount  of  viscous  response,  and  relatively  little  plastic 
response.  As  the  stress  level  increases,  plastic  response  becomes  more 
and  more  significant. 

The  second  example,  sea-ice,  illustrates  a highly  viscous  material 
with  pronounced  plasticity  and  relatively  small  elastic  response.  The 
experimental  data  were  generated  by  Vaudrey  [28]  and  are  shown  as  dis- 
crete data  points  in  Figure  19.  The  data  shown  are  for  sea-ice  at  -27°C, 
with  ice  crystals  orthogonal  to  the  direction  of  loading  at  three  com- 
pressive stress  levels,  c = 70,  175,  and  350  psi. 

A viscoelastic-plastic  model  fit  gives  the  expression  for  axial 
strain  as: 
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It  is  assumed  c.  = 250  psi,  and  the  hardening  function,  E (e  ) * o.  / 

2 L P P L 

(1  + c ) , which  approaches  zero  as  t becomes  large. 

P P 

The  above  model  expression  is  plotted  as  solid  lines  in  Figure  19, 
showing  excellent  agreement  with  the  data  points.  In  particular,  the 
viscoelastic-plastic  model  is  able  to  replicate  unrestrained  plastic 
flow  that  occurs  at  the  high  stress  level. 


o = 350  psi  (instantaneous  plastic  flow) 
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Figure  19.  Comparison  of  viscoelastic-plastic  model 
with  experimental  data  for  sea-ice. 
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SUMMARY  AND  RECOMMENDATIONS 

This  report  offers  a textbook  style  development  of  viscoplasticity. 
For  completeness,  a detailed  review  of  viscoelasticity  and  plasticity  is 
also  included.  The  main  text  focuses  on  the  so-called  "viscoelastic- 
plastic  model."  However,  this  model  is  just  one  example  of  a general 
family  of  combo-viscoplastic  models  presented  in  the  appendixes. 

Each  model  is  initially  presented  with  visual  one-dimensional 
concepts  and  then  generalized  into  multidimensional  stress  space.  The 
net  result  is  an  incremental  constitutive  model  suitable  for  numerical 
computation.  To  this  end,  finite  element  algorithms  are  proposed  for 
each  model. 

The  viscoelastic-plastic  model  is  shown  to  capture  the  observed 
creep  responses  of  such  elusive  materials  as  plexiglass  and  sea-ice. 
Also,  it  is  apparent  other  materials,  such  as  soil,  plastic,  epoxy, 
concrete,  etc.,  can  be  suitably  approximated  with  this  model  within 
certain  ranges. 

A primary  advantage  of  the  viscoelastic-plastic  model  is  the  ease 
of  parameter  identification.  This  is  because  it  can  be  divided  into 
separate  viscoelastic  and  plastic  identification  problems.  The  primary 
limitation  is  that  the  strain  rate  is  independent  of  plastic  strain  for 
a constant  load.  This  limitation  can  be  removed  by  employing  higher 
order  combo-viscoplastic  models. 

Immediate  future  efforts  should  be  directed  on  two  fronts.  First, 
the  algorittims  developed  herein  should  be  incorporated  into  a finite 
element  program  and  tested  on  selected  boundary  value  problems.  Second, 
and  concurrently,  higher  order  combo-viscoplastic  models  should  be 
investigated  and  compared  with  experimental  test  data.  Here  t ho  objec- 
tive is  to  identify  a single  "grandfather"  combo-vi. scop  las t ic  model  that 
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contains  all  the  desired  response  characteristics  and  is  capable  of 
degenerating  to  lower  order  models. 

The  net  result  of  these  efforts  would  be  the  beginnings  of  an 
analytical  tool  (finite  element  program)  capable  of  rationally  dealing 
with  structural  engineering  problems  encountered  by  the  Navy. 

Finally,  the  long-range  objective  should  include  a material 
identification  study.  That  is,  for  each  material  type,  subtypes  would 
be  grouped  by  ranges  of  state  variables,  such  as  temperature,  moisture, 
impurities,  etc.  Then,  visceplastic  parameters  (i.e.,  parameters  of  the 
"grandfather"  model)  could  be  determined  by  identification  procedures 
discussed  herein.  This  library  of  data  could  be  incorporated  into  the 
finite  element  program  so  that  model  parameters  would  be  automatically 
retrieved  (or  interpolated)  by  specifying  the  material  type  and  state 
variables.  Alternatively,  experimental  data  (e.g.,  triaxial  test  data) 
could  be  input  into  the  program  and  the  parameters  determined  directly. 

The  benefits  of  this  plan  are  enormous,  both  technically  and 
economically.  On  the  technical  side  the  plan  provides  a rational  and 
uniform  approach  for  dealing  with  structural  materials,  such  as  concrete, 
plastics,  soils,  sea-ice,  epoxies,  and  a host  of  other  time-dependent 
materials  encountered  day-to-day  by  Navy  engineers.  Economically,  under 
current  procedures,  it  is  not  uncommon  for  an  engineer  to  spend  more 
than  one-half  of  his  time  groping  with  material  models  and  their  para'n- 
eters.  The  proposed  plan  would  relieve  the  engineer  of  this  burden  and 
allow  him  to  get  on  with  solving  the  problem. 
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Appendix  A 


PLASTICITY  MODELS  AND  CONCEPTS 


In  this  appendix,  plasticity  concepts  are  elaborated  for  common 
yield  conditions  and  flow  rules  with  emphasis  on  methods  of  identifying 
plasticity  parameters  from  experimental  data. 

The  first  and  foremost  requirement  to  establish  a plasticity  model 
is  the  selection  of  a yield  condition.  UBing  the  universal  form  estab- 
lished in  Equation  2-34,  the  general  yield  condition  is  written  as: 


F - f(o  - o ) - k*  (A-l) 

~ P 

such  that  F 
F 

and  f(o  - o ) 

**  P 

k* 


= 0 implies  plastic  response 

< 0 implies  nonplastic  response 

= loading  function 

= yield  parameter  (universal  weighting) 


For  a particular  current  value  of  the  weighted  plastic-tracking 
stress,  5p,  and  of  the  yield  parameter,  k* , the  yield  surface  is  defined 
by  all  stress  states,  o,  such  that  F = 0.  The  unit  normal  of  the  yield 
surface  in  six-dimensional  stress  space  is  given  by: 
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where 


3 f (o  - o ) 


The  unit  normal,  n,  is  an  important  vector,  because  it  controls  the 
T /N  ^ 

magnitude  (Ao  • n)  and  direction  of  plastic  strain.  Observe  from 
Equation  A-2  that  n is  determined  directly  from  the  gradient  of  the 

*v 

loading  function  so  that  the  choice  of  loading  function  dictates  the 
general  nature  of  plastic  flow. 

Loading  functions  are  generally  written  as  functions  of  the  stress 
invarients: 


f(o  - 5 ) 

**  <s*  p 


f l 0 (o  - o ),  j'(o  - a ),  J ' (n  - a ) 

IT)  **  P Z <v  J • 


(A-4) 


where  o is  the  spherical  stress  invariant  and  J'  and  J'  are  the  second 
m r 2 3 

and  third  deviatoric  stress  invarients.  Using  the  chain  rule,  the 
gradient  vector  is  given  by: 


- 3f  ^ 3f  Sf 

9o  ~1  9 J L ~2  9J  ' ~3 

m 2 3 


(A-5) 
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The  vectors  a,,  a , and  a are  independent  of  the  loading  function 
and  are  given  by: 


3o 
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~1  do 


(A-6) 
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The  foregoing  has  illustrated  general  forms  for  yield  conditions 
with  arbitrary  loading  functions.  In  the  next  section  these  concepts 
will  be  illustrated  on  a particular  form. 


DRUCKER-PRAGER  YIELD  CONDITION 

A particularly  useful  and  versatile  yield  condition  is  the  so-called 

Drucker-Frager  model.  For  clarity  and  conciseness,  the  model  will  be 

developed  assuming  isotropic  hardening  (i.e.,  o = 0);  however,  the 

P 

inclusion  of  universal  hardening  is  a straightforward  extension. 

The  Drucker-Prager  yield  condition  is: 


F = 


(A-9) 
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which  implies  the  loading  function  is: 


f = “ °m  + V7] 


(A-10) 


The  invar ients  and  are  measures  of  hydrostatic  stress  and 
shear  stress,  respectively,  and  are  defined  as: 


o 

m 


1 
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(°1  + o2  + o ) 


(A-ll) 


J2  ■ T l°l  ' °»)2  + (°2  - °m>2  + (03  - V2' 


. 2 2 2 
+ T12  + T23  + T31 


(A-12) 


where  o^,  and  are  normal  stresses,  and  T23’  anc*  T31  are 

shear  stresses. 

The  parameters  of  the  yield  condition  are  a and  k*,  which  are 
selected  to  establish  the  shape  of  the  yield  surface.  For  example,  if 
a = 0,  the  yield  surface  becomes  a cylinder  of  revolution  about  the 
hydrostat  in  principal  stress  space  with  radius  k*,  or  equivalently,  it 
may  be  viewed  as  a straight  line  on  a o , & graph  as  shown  in  Figure 
A-la.  This  type  of  loading  function  is  associated  with  the  Von  Mises 
yield  criterion  and  is  a good  representation  f or  ductile  metals. 

For  a more  general  case,  -i  is  specified  non  zero,  which  provides 
the  conical  surface  shown  in  Figure  A-lb.  The  conical  surface  is 
associated  with  the  standard  Drucker-Prager  yield  criterion,  which  has 
applications  in  concrete,  soil,  rock,  etc. 


°1  = °2  = °3 


(b)  Drucker-Prager  yield  criterion. 


The  yield  condition  can  also  be  applied  in  piecewise  continuous 

fashion  to  better  approximate  actual  yield  surfaces  as  indicated  in 

Figure  A-lc.  To  apply  this  model,  one  loading  function  is  operable  for 

values  of  c < o , and  a second  loading  function  is  operable  for  a > c . 

mo  mo 

Additional  functions  can  be  added  if  desired.  This  is  termed  "modified 
Drucker-Prager"  model. 

The  Von  Mises  yield  criterion  only  requires  one  parameter,  k*,  and 
therefore,  only  one  material  test  to  initially  characterize  the  yield 
surface.  The  standard  ^rucker-Prager  criterion  requires  determination 
of  a and  k* ; thus,  a r ’ two  material  tests  are  required.  Table  A-l 
identifies  the  Vo.  '■  >.<;  ?■  meter  for  a tensile  test  and  the  Drucker- 

Prager  parameters  >■  ' ^ . iXial  strength  test. 

The  description  . and  k*  for  Drucker-Prager  are  in  terms  of  the 
cohesion  stress,  C,  and  internal  angle  of  friction,  6,  as  determined 
from  at  least  two  triaxial  tests.  To  illustrate  the  derivation  of  a and 
k*,  consider  two  triaxial  tests  under  different  confining  pressures  and 
loaded  axially  to  failure  as  suggested  in  Figure  A-2.  Drawing  Mohr 
circles  for  both  stress  states  at  failure  permits  defining  a straight 
line  failure  surface  with  constants  C and  6 as  illustrated  in  Figure 
A-3.  The  failure  surface  prescribes  what  combination  of  shear  stress 
and  normal  stress  will  cause  failure,  i.e.  , 


t = C + o tan  t 


(A- 131 


where  : = shear  stress  at  tailure  (positive) 

normal  stress  at  failure  (compression  positive) 


All  states  of  stress  (i,  •••■)  satisfying  Equation  A-13  are  failure 
states,  and  the  equation  is  known  as  a Coulomb  failure  criterion.  The 
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rotated  stresses,  *.  and  a,  are  related  to  the  principal  stresses,  and 
o ^ , through  the  Mohr  circle  geometry  as: 


cos  c 


(A-14) 


o 


2 


sin  *3 


(A- 15) 


Accordingly,  Equation  A-13  may  be  written  as: 


_ / 1 + sin  0 \ , 

l 1 - sin  u j - 


2 C cos 

1 - sin 


= 0 


( A — 1 6 ) 


Now  the  objective  is  to  put  the  general  yield  condition,  a c + 

m 

- k*  = 0,  into  the  form  of  Equation  A-16  in  order  to  equate  co.  f i- 
cients  and,  thereby,  determine  k*  and  a in  terms  of  C and  0.  For  a 
triaxial  test,  we  have  = -(r ^ + 2 r^)/3  and  - r^)/^T 

(assuming  compression  positive  and  ^ > • ^) , so  that  the  general  yield 
condition  can  be  written  as: 


(VI  + 2 ct) 
(VI  - 0 


3 k* 


vi  - , 


= 0 


( A - 1 7 ) 


Equating  the  coefficients  of  ; ^ and  the  constant 
Equations  A-16  and  A-17  gives  ttie  values  for  k*  and  -i 


terms  between 
in  Table  A-l. 
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It  should  be  appreciated  that  the  determination  of  a and  k*  can  be 
achieved  in  a direct  manner  without  the  need  of  computing  C and  6. 

That  is.  Equation  A-17  can  be  used  directly  by  inserting  the  failure 
stresses,  and  for  the  two  test  specimens,  thereby  providing  two 
equations  for  the  two  unknowns  k*  and  a. 


DRUCKER-PRAGER  ASSOCIATIVE  FLOW  RULE  AND  HARDENING  FUNCTION 

In  general,  an  associative  rule  (e.g.,  Equation  2-26)  can  be 
written  as: 


Ac 

-P 


H’ 


N Ao 


(A-18) 


where  N 


A aT 


n n 


In  the  above,  Ac  is  the  plastic  strain  increment,  Aa  is  the  total 

~p  ~ 

stress  increment,  H'  is  the  hardening  function,  and  N is  dependent  on 
the  unit  normal  of  the  yield  surface,  n. 

Consider  the  Dr ucker-Prager  loading  function  (Equation  A-10)  for 
states  of  stress  characteristic  of  common  tests  (such  as,  triaxial  or 
confined  compression);  i.e.. 


independent  stress 


'J  2 U3: 

~ 1 2 = ”'13 


independent  stress 


= i 


23 


US 


0 


The  yield  surface  unit  normal  for  this  stress  state  is 


(A-20) 


where  6 


A 


sgn 


2 a -!-  2VT  sgn 
2 a - yT  sgn 

2 a - VT  sgn 

\IV(2  -J  + 3) 

sign  of  (o^  - Oj) 


Assuming  the  material  is  initially  isotropic,  the  two  independent 
stresses  (say  and  c^)  produce  a corresponding  strain  state  as  follows: 

independent  strain 

l 2 = 0^:  independent  strain 

V12  = Y23  = *13  = 0 

With  the  above  symmetry,  the  plastic  flow  rule  (Equation  A-18)  can 
be  reduced  to  the  following: 


plastic 


total 


(A-21) 
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The  corresponding  elastic  strains  for  this  stress  3tate  are: 


A£1 

Ae, 


elastic 


1 -2  v 

-v  (1  - 


“’ll 

rn 

D 

< 

(A-22) 


total 


Therefore,  the  total  stress-strain  relationship  is  given  oy  the  sum  of 
elastic  and  plastic  strain  as: 


ACj 

Ae, 


/ A2  62  . iW  2 A2  g 2 v\ 

\ H'  E/V  H'  E 7 


a2  e 

H' 


v 

E 


A0] 

Ac, 


(A-23) 


Equation  A-23  is  applicable  for  any  test  data  conforming  to  the  symmetry 
conditions:  Ae 2 * Ae^i  Ao2  = Ao^,  and  all  shear  stress  and  strains  are 
zero. 

For  a triaxial  test,  Ao^  and  A are  specified,  and  Ae^  is  the 
measured  axial  strain  increment  (elastic  and  plastic)  . Assuming  the 
elastic  properties  E and  v are  determined  by  standard  procedures,  the 
hardening  function  H'(e^)  can  be  determined  from  Equation  A-23  as: 


H ' ( c ) 
P 


A B ( 8 Aax  + 2 Ao2> 

Ao  2 v Jo 

Aci  - — + — ¥— 


(A-24) 
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As  another  example,  a one-dlmenslonal  strain  test  specifies  Ao^  and 
measures  Ae^  with  the  constraint  Ae^  ■ 0.  For  this  case,  the  hardening 
function  is  given  by: 


H'(ep) 


A2[2  + 82(1  - v)  + 4 8 v]  - (2  A2  E) 

1 2 Ael 

4(2  v2  + v - !)  (1  - v) 


(A-25) 


It  is  interesting  to  note  that  for  the  one-dimensional  strain  case, 
a perfectly  plastic  material,  H^e^)  » 0,  does  not  result  in  unrestrained 
plastic  flow.  This  is  because  the  lateral  stress,  o^,  adjusts  to  keep 
the  stress  state  on  the  yield  surface. 


DRUCKER'S  POSTULATE  AND  FLOW  RULES 


Drucker's  postulate  may  be  stated  as:  "Given  a stress  state  in 
equilibrium  on  the  yield  surface  and  an  external  load  cycle  slowly 
applied  and  removed,  then  the  external  agency  does  positive  work  during 
loading,  and  the  net  work  over  the  cycle  is  non-negative." 

A material  conforming  to  this  postulate  is  termed  stable.  Also,  a 
mathematical  flow  law  ensuring  this  condition  is  admissible  (i.e. , will 
not  create  energy).  In  equation  form,  Drucker's  postulate  is: 


do  (de  + de  ) > 


(loading) 


do^(dc  + dc  ) 
» „p  -e 


doT(dr  ) > 0 


(full  cycle) 
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"TP 


i 

i 


The  £irst  relationship  signifies  both  elastic  and  plastic  work 
increments  are  positive  during  loading.  The  second  relationship  implies 
unloading  work  is  all  elastic  so  that  for  a full  cycle  we  must  have: 


do^  dt  > 0 

~P  - 


(A-26) 


where  m is  the  assumed  direction  of  plastic  straining.  If  m = n,  the 

flow  rule  is  said  to  be  associative;  otherwise,  it  is  nonassociative. 

T A 

By  virtue  of  the  vield  criterion,  the  scalar  product  do  . n is 

always  non-negative  during  plastic  deformation.  Therefore,  Equation 

A-27  is  always  satisfied  for  an  associative  flow  rule,  since  H'  is  also 

non-negative  for  stable  materials. 

In  order  for  a nonassor iat ive  flow  rule  to  satisfy  the  stability 

T ^ 

criterion,  we  must  have  do  • m > 0 for  any  do  that  causes  plastic  flow. 
Further  implications  of  nonassoc iative  rules  are  beyond  the  scope  of  this 
appendix  and  require  additional  research. 
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Appendix  B 


IDENTIFICATION  OF  VISCOELASTIC  RELAXATION  AND  CREEP  FUNCTIONS 


The  intent  of  this  appendix  is  to  provide  guidelines  for  converting 
viscoelastic  data  into  creep  and/or  relaxation  functions  of  an  exponen- 
tial series  form.  Specifically,  the  following  three  basic  problems  will 
be  addressed: 

1.  Obtaining  a creep  function  from  creep  data,  or  obtaining 
a relaxation  function  from  relaxation  data. 

2.  Inverting  a creep  function  to  a relaxation  function  and 
vice  versa. 

3.  Transforming  relaxation  functions  from  one  type  to  another, 
such  as  transforming  a Young's  Modulus  function  into  bulk 
and/or  shear  functions. 


I.  CREEP  AND  RELAXATION  FUNCTIONS  FROM  VISCOELASTIC  DATA 

The  objective  is  to  curve-fit  (by  least-square  error)  an  exponen- 
tial series  with  available  viscoelastic  data  to  achieve  either  a creep 
or  relaxtion  function.  The  exponential  series  has  the  general  form: 


N -X  t 

Y ( t , 1 ) = A + B t + £ C,  (1  - e ) (B-l ) 

k=l 
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where  A,  B.  C^,  Cj,  ...  C^,  X^,  Xj,  ...  X^  are  unknown  parameters  to  be 
determined.  In  this  form,  the  term  A is  the  initial  elastic  response,  B 
is  the  rate  of  linear  response,  is  the  long-time  viscous  response, 
and  X^  is  the  inverse  of  retardation  or  relaxation  time.  For  conveni- 
ence, these  parameters  are  denoted  by  <fr,  i.e. : 


4> 


2 ’ 


M 


A , B , ^ , Cj, 


where  M = 2(N  + 1) 


The  form  of  Equation  B-l  is  applicable  for  either  a creep  or 
relaxation  function  with  the  necessary  provisions  of  Y(t,t>)  > 0 and 
> 0 to  insure  decay  of  each  exponential  term.  In  the  case  of  a 
relaxation  function,  Y(t,$)  must  be  a monotonically  decreasing  function 
so  that  the  linear  parameter  B = 0.  However,  a creep  function  mono- 
tonically increases;  thus,  the  parameter  B need  not  be  specified  zero. 
Figure  B-l  illustrates  example  forms  of  Equation  B-l  for  creep  and 
relaxation  functions. 

The  known  viscoelastic  data  represent  either  creep  data  or  relaxa- 
tion data.  Here,  creep  data  are  defined  as  a recorded  strain  history 
corresponding  to  a material  specimen  subjected  to  a constant  unit  of 
stress,  such  as  axial  stress,  shear  stress,  or  hydrostatic  stress. 
Conversely,  relaxation  data  are  defined  us  p.  recorded  stress  history 
corresponding  to  a material  specimen  subjected  to  a constant  unit  of 
strain  deformation  (c.g.,  axial,  shear,  volumetric,  etc.) 

In  either  case,  it  is  assumed  the  data  are  reduced  to  a set  of 
discrete  values  (Y  ,t^),  i = 1,  2,  ...  n,  where  Y^  is  the  creep/ 
relaxation  data  point  at  time  t^. 
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These  discrete  data  may  be  written  as  a piecewise  linear  function 
of  time  as  follows: 


Y(t)  = V±  + t,  for  t±  < t < t1+1  (B-3) 


where  W , 


(Yi  ' Yi  V 


(Vi+1  - *,>'<‘1+1  - t,> 


In  the  above,  Y(t)  is  a piecewise  function  connecting  successive  data 
points,  as  illustrated  in  Figure  B-2  for  the  case  of  creep  data.  Higher 
order  interpolation  functions  could  be  used  if  desired. 

To  fit  the  function  Y (t,<J>)  (Equation  B-l)  to  the  viscoelastic  data 
Y(t)  (Equation  B-3),  the  least-square-error  method  with  Newton-Raphson 
solution  procedure  is  used  as  follows. 

First  define  error: 


error(t,4>)  = Y(t,$)  - Y(t) 


(B-4) 


Squaring  the  error  and  summing  it  by  integrating  over  the  time  of 
interest,  0 _<  t <_  tn,  gives  the  net  square  error: 


e(4>) 


(Y(t,<f)  - Y ( t ) ] 2 dt 


(B-3) 
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Minimizing  this  net  error  with  respect  to  each  unknown  parameter, 
<(>  = (A,  B,  C,  ...  C^,  , X2,  ...  X^)  gives  M equations: 


3c(<fr) 


= 2 


/ 


[Y(t,(J>)  - Y ( t)  ] 


m 


3Y(t,») 

3<> 


dt  = 0 


(B-6) 


where  m = 1,  2,  . . . M.  For  brevity,  the  above  equations  are  denoted  by 
an  M-ditnensional  vector  equation  as: 


£ (<J>0)  = 0 


(B-7) 


where  = |a°,  B°»  ck°’  ] represents  the  solution  vector  satisfying 
Equation  B-7.  Since  the  equations  are  nonlinear  with  respect  to  the 
parameters  X^»  an  iterative  solution  technique  is  required.  In  the 
following,  a Newton-Raphson  iterative  solution  procedure  is  presented. 

The  Newton-Raphson  procedure  is  based  on  a Taylor  series  expansion 
of  Equation  B-7,  where  $°  is  written  as: 


(B-8) 


s o 

such  that  i is  any  estimate  of  $ , and  At  corrects  the 
satisfy  Equation  B-8.  The  objective  is  to  successively 

s 

estimate,  of  i so  that  A f ->  0.  To  this  end,  the  Taylor 
of  Equation  B-7  can  be  written  as: 


estimate  to 
obtain  a better 
series  expansion 
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0 


f(4>S  + At)  = 0 = f(4>3)  + F (<J>S ) • At  + higher  order  terms  (B-9) 


where  F(t  ) = 3(f)/3t 


S 8 

F(t  ) is  an  MxM-dimensional  Jacobian  matrix  evaluated  at  t , and 
the  higher  order  terms  are  second  order  and  above  products  of  A<t>  (i.e., 

A# 

squares,  cubes,  etc.). 

If  At  is  small  so  that  higher  order  terms  are  negligible,  then  At 

a# 

may  be  approximately  determined  as: 


At  - -F-1(tS)  • f(t9) 
— ■ •.  - ~ 


(B-10) 


With  this  approximation  for  At  a better  estimate  for 


Q 

t is: 

A# 


S+l  8 

<p  * <t>  + A<|> 

M W 


(B-ll) 


S+1 

Returning  to  Equation  B-10,  the  process  is  repeated  utilizing  t in 

3 

place  of  t • After  a sufficient  number  of  iterations,  we  have  At  0 

8+1  "**0  ** 
and  t -*■  t » providing  the  desired  solution. 

A*  A* 

The  major  effort  in  the  above  algorithm  is  establishing  the  vector 
s s 

f(t  ) and  the  Jacobian  F(t  ) for  each  iteration.  These  somewhat  labo- 

a m a* 

rious  developments  are  given  next. 

For  specificity,  the  M equations  represented  by  f(t)  = 0 and 

^ A* 

defined  by  Equation  E-6  are  distinguished  by  the  parameter  groups  A,  B, 
Cj , and  Xj  as  follows: 
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[Y  (t,  (J>)  - Y(t)  ] dt 


(B-12  ) 


£ 


B 


3e(&) 

3B 


A 

J [Y(t,4>)  - Y(  t)  ] 
'A 


t dt 


(B— 1 3) 


fr  * ^§7^-  = f ~ V(t)](l  - e j ) dt 

cj  3Cj  0 

j = 1,  2,  ...  N (B-14) 


ae(t) 

3XJ 


- -Xifc 

J (Y(t , 4>)  ~ Y(t)  ] C t e 3 dt 
0 3 

j = 1,  2,  . . . N (B-15) 


where  the  common  factor  2 has  been  divided  out  of  all  M equations. 
Performing  the  indicated  time  integration  for  Equations  B-12  through 
B-15  gives  the  results  presented  in  Table  B-2,  where  the  functions  P,  Q, 
R,  and  S are  defined  in  Table  B-l. 

Next,  expressing  the  components  of  the  Jacobian  matrix  F in  a 

* 

similar  notation,  we  have: 


F 


AA 


(B-16) 
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Table  B-l.  Function  Definitions  and  Values  for  Some  Exponential  Integrations 


A &&&  «.  ifllltt'iiMft L m • r iliiiiiiiiiNItifliiiiii.i ■■ . ii a i.JilLi., 'llUiiJKkii-  A'rJ- 


The  integrated  expressions  for  Equations  B-16  through  B-25  are  given  in 

Table  B-3.  Note  F is  always  symmetric  by  virtue  of  mixed  partial  deriva- 

" ? 

tives,  i.e.,  F$  = 3 e/3^  3<’j  • 

To  summarize,  the  following  steps  comprise  the  curve-fitting 
algorithm: 

1.  Make  initial  estimate  of  parameters: 


=_  Q 3 

P 2.  Compute  error  functions  f($  ) and  Jacobian  matrix  F(4>  ) 

| (Tables  B-2  and  B-3) . 

f 1 

? 3.  Solve  for  estimate  correction  A<J>,  i.e.. 


WrT'rwfVT!*'^  WTr'F'P'H' I"i‘l''",",*|:",n":')"1  •"l"1"1  ' •"""  '■  * 7 ~~l  y?  rs  W1|*PW?' 1 


S+l  S I | 

4.  Update  estimate,  <f  =4)  + A4>;  if  | A<J> j is  near  zero,  solution 

8^*1 

is  complete.  Otherwise,  return  to  step  2 using  4>  to  replace 


For  the  implementation  of  this  algorithm  it  is  prudent  to  provide 
the  capability  of  specifying  any  combination  of  the  parameters  at  pre- 
determined values.  This  is  useful  if  the  algorithm  cannot  converge  on  a 
solution  when  all  parameters  are  unspecified.  In  particular,  the  con- 
stant term  A should  be  specified  as  the  initial  elastic  response.  This 

-X  t 

Is  the  reason  the  exponential  terms  were  written  as  C(1  - e ),  instead 
— Xt 

of  simply  Ce  as  presented  in  the  main  text.  Also,  if  all  X^  are 
specified,  then  a linear  solution  is  achieved  (no  iterations). 


II.  INVERTING  CREEP  AND  RELAXATION  FUNCTIONS 

Once  a creep  or  relaxation  is  determined,  its  inverse  can  be 
determined  by  the  procedure  outlined  in  this  section.  It  Is  assumed  the 
creep  or  relaxation  functions  are  characterized  by  exponential  series  as 
i ollows : 


M -B  . c 

J(t)  =J  +JLt+y'j.eJ  (B-26) 

3 b j-1  J 


Y(t) 


Y + 
a 


N -cia  t 

Iv 

i = l 


(B-27) 


In  the  abcve,  J(t)  represents  a general  creep  function,  and  Y(t) 
represents  a general  relaxation  function.  These  forms  intentionally 
exclude  materials  that  do  not  exhibit  an  initial  elastic  response. 
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The  objective  can  now  be  stated  as:  given  the  parameters  of  the 

creep  function  (J  , J,,  J, , ...  Ju,  3.,  ...  8W)  , what  are  the  corre- 
a d t Mi  M 

sponding  parameters  of  relaxation  function  (Y  , Y^,  ...  Y^,  a^,  ...  aN)? 
Or  conversely,  given  the  relaxation  parameters,  what  are  creep  param- 
eters? 

As  illustrated  in  Figure  B-l,  the  behavior  of  a material  is 
different  depending  on  whether  the  material  is  solid-like  or  fluid-like. 
Accordingly,  the  creep  and  relaxation  functions  may  be  subdivided  into 
two  types  as  follows: 

1.  Solid-like  creep/relaxation  function  forms  (M  * N,  = 0) : 


J(t)  - J + £ J. 


(B-28) 


Y(t)  - Y„ 


+ E Yi e 

i-i 


(B-29) 


2.  Fluid-like  creep/relaxation  function  forms  (M  ■ N - 1): 


N-l  -8 . t 


J + J.  t + D J.  e ^ 
3 b j - 1 J 


(B-30) 


N -a , t 

ZV 

i-1 


(B-31) 


where  the  following  restrictions  on  the  parameters  are  implied: 


J 


V 


> o 


j 


t 


< 0 


J + 
a 


M 

L ^ 


j-i 


> 

i 


0 


Note  the  solid-like  creep  and  relaxation  functions  have  identical  forms: 
each  contain  a constant  term  and  decaying  exponential  series.  However, 
the  fluid-like  creep  function  contains  the  linear  term  J^t  which  permits 
continued  creep  (or  flow)  for  all  time.  The  corresponding  fluid-like 
relaxation  function  contains  no  constant  term  or  linear  term,  but  rather 
an  additional  exponential  term.  This  permits  the  stress  to  relax  to 
zero  under  a fixed  strain. 

To  determine  the  respective  inverse  functions  for  either  the  fluid- 
like or  solid-like  forms,  the  fundamental  convolution  identity  between 
creep  and  relaxation  functions  is  exploited: 


J*  dY 


- 1 


or,  in  expanded  form: 
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1 


(B-32) 


J(t)  Y{0)  + 


J(t  - T)  -dT 


For  the  solid-like  case  each  function  contains  2N+1  parameters. 

Upon  inserting  Equations  B-28  and  B-29  into  the  above  identity  and 

-a^t 

integrating,  the  results  may  be  grouped  into  time  functions  of  e , 
-Sit 

e , and  a constant.  The  groups  of  terms  associated  with  each  expo- 
nential time  function  must  be  zero , and  the.  constant  tennis  must  equal  1 
in  order  to  satisfy  Equation  B-32  for  all  time.  This  gives  2N+1  rela- 
tionships between  the  constants  as  follows: 


J Y =1 
a a 


(B-33) 


Ja  £ (6k  - V 


«i  L J,  TT  (B 

j=l  J k-1 


i « 1.  2, 


(B-34) 


ir  (B 

k=l  J 


1,  \ Yi  7T 

i=l  k=l 


j = 1,  2, 


(B-33) 


If  the  creep  function  parameters  are  known,  the  relaxation  parameters 


V V •• 


uM  are  determined  by  solving  for  the  roots  of  the  N-degree 


polynomial  of  u in  Equation  B-34.  Next,  the  relaxation  coefficients  Yj , 
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Y^,  . . . Y^  are  determined  by  the  linear  solution  of  Equation  B-35,  with 

Y = 1/J  . On  the  othe-  hand,  if  relaxation  parameters  are  known,  the 
a a 

creep  parameters  8^,  82.  •••  6^  ate  determined  from  the  roots  of  Equa- 
tion B-35  first,  then  the  creep  coefficients  are  determined  from  Equation 
B-34,  with  J = 1/Y  . Table  B-4  provides  solutions  for  N = 1 and  2. 

In  a similar  manner,  parameter  relationships  for  fluid-like  mate- 
rials can  be  developed  by  inserting  Equations  B-30  and  B-31  into  the 
convolution  identity.  For  convenience,  the  following  expression  is 
used . 

f Y*  dJ  = 1 

*f) 


For  this  case,  there  are  2N  parameters  per  function  along  with  2N 
relationships  between  the  creep  and  relaxation  functions  given  as: 


'.tall"'.  WLMiij.1,1':.. 


Relaxation  Function  Inversions  for  Solid-Like  Materials 


N N 

E Yi  ^ 

i=i 


(<».  - 


k=i 

k*i 


V 


= 0 


j = 1,  2, 


N-l 


(B-38) 


If  creep  parameters  are  known,  the  relaxation  parameters  cx^  a,,,  ... 
are  determined  from  solving  the  roots  of  che  polynomial  given  by  Equa- 
tion B-37.  Thereafter,  Y^  Y v ...  Y^  are  determined  by  the  linear 
solution  of  the  N-l  equations  of  Equation  set  B-38  together  with  Equation 
B-36.  On  the  other  hand,  if  relaxation  parameters  are  given,  the  N-l 
roots  of  Equation  B-38  give  8^  Br  •••  BN_r  and  the  creep  coefficients 
J J , J , Jj,  •••  Jjq  ^ are  determined  by  the  N equations  in  Equation 
set  B-37  together  with  Equation  B-36.  Table  B-5  shows  solutions  for 
N = 1 and  2. 


III.  RELAXATION  FUNCTION  INTERRELATIONSHIPS 

Thus  far,  relaxation  and  creep  functions  have  been  denoted  in  a 
general  manner  by  the  functions  Y(t)  and  J(t),  respectively.  However, 
in  practice,  Y(t)  represents  a particular  kind  of  relaxation  function, 
such  as  shear  modulus,  bulk  modulus,  Young's  modulus,  or  confined 
modulus,  etc.  Similarly,  J(t)  represents  the  corresponding  inverse 

function. 

If  experimental  data  are  available  for  one  particular  type  of 
function,  it  may  be  necessary  to  convert  the  data  to  another  type  of 
function.  For  example,  suppose  a Young's  modulus  relaxation  function 
E(t)  has  been  determined  from  a test  coupon  along  with  an  observed  time 
function  for  Poisson’s  ratio  v(t).  If  it  is  desired  to  determine  the 
equivalent  bulk  and  shear  relaxation  functions,  then  the  material  data 
must  be  transformed. 


liJ 


if 


where  D , D,  , and  D are  functions  of  time,  and  the  svmbol  * is  the 

convolution  operator.  It  is  evident  that  Equation  B-39  has  the  same 

form  as  an  elastic  constitutive  matrix,  and,  like  dn  elastic  matrix,  the 

terms  D^,  D^,  and  are  composed  of,  at  most,  two  independent  material 

properties,  which  in  this  case  are  functions  of  time.  Possible  time 

function  pairs  are:  Young's  modulus  and  Poisson's  ratio,  bulk  modulus 

and  shear  modulus,  and  confined  modulus  and  lateral  stress  coefficient. 

Clearly,  it  makes  no  difference  what  pair  or  combination  of  material 

time  functions  is  used  as  long  as  they  combine  to  produce  the  same  time 

functions  for  D (t),  D,(t),  and  D (t).  Therefore,  the  various  visco- 
a b c 

elastic  moduli  (functions)  have  the  same  interrelationships  as  the 
corresponding  elastic  moduli. 

Table  B-6  provides  some  useful  interrelationships  to  aid  in  trans- 
forming from  one  type  of  function  to  another. 
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Table  B-6.  Interrelationships  Among  Elastic  and 
Viscoelastic  Material  Functions 


Young's  Modulus 

Bulk  Modulus 

Confined  Modulus 

Material 

and 

and 

and 

Function* 

Poisson's  Ratio, 

Shear  Modulus, 

Late.al  Coefficient, 

E,  v 

K,  G 

M , X 
s ' 0 

Young's 

Modulus, 

E - 

1 

E 

9 G K 

(1  + 2K  ) (1  - K ) 

M o O 

3K  + C 

MS  1 + K 

0 

Poisson' s 
Ratio, 

V “ 

3K  - 2G 

K 

o 

V 

2 (3K  + G) 

1 + K 

o 

Bulk 

M 

Modulus , 

K * 

£ 

K 

-f  11  + 2V 

3(1  - 2v) 

Shear 

M 

Modulus, 

G - 

£ 

•f (1  - V 

2(1  + v) 

0 1 

Confined 
Modulus , 

M - 

s 

E(1  -v) 

K +i-G 

M 

s 

(1  + v)(l  - 2v) 

Lateral 

Coefficient, 

Ko- 

V 

3K  - 2G 

K 

o 

1 - V 

3K  + 4G 

*A11  parameters  may  be  functions  of  time. 
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To  illustrate  the  procedure  for  transforming  viscoelastic  functions, 
consider  the  example  of  determining  bulk  and  shear  relaxation  functions 
given  test  data  from  a simple  creep  test.  That  is,  a test  coupon  of  the 
material  is  axially  loaded  with  a constant  stress,  and  axial  and  radial 
strain  are  recorded.  Then  the  following  steps  are  required: 

1.  Normalize  experimental  data.  Divide  axial  strain  by  applied 
stress  to  get  creep  data  (this  corresponds  to  inverse  of  Young's  modulus 
relaxation  function).  Also,  divide  radial  strain  data  by  axial  strain 
data  to  get  Poisson's  ratio  function,  v(t). 

2.  Determine  creep  function.  Using  the  creep  data,  employ  the 
methods  of  Part  I of  this  Appendix  to  determine  the  best-fitting  creep 
function. 

3.  Determine  Young's  modulus  relaxation  function.  Invert  the 
above  creep  function  to  get  the  corresponding  Young's  modulus  relaxation 
function,  E(t) . Use  the  techniques  in  Part  II  of  this  Appendix  for 
inverting. 

4.  Transform  to  bulk  modulus  relaxation  function.  From  Table  B-6, 
the  bulk  modulus  is  given  as  K(t)  = 1/3  E(t)/[1  - v(t)].  If  the  experi- 
mental data  show  v(t)  is  constant  (or  near  constant),  then  K(t)  has  the 
same  form  aa  E(t)  and  only  differs  by  the  scalar  divisor,  3(1  - v) . 

If,  on  the  other  hand,  v(t)  is  not  constant,  K(t)  can  be  determined  by 
the.  best-fit  procedure  in  Part  I of  this  Appendix  where  the  data  points 
are  Kft^  = 1/3  E(t)1/[1  - v ( t ±)  ] . 

5.  Transform  to  shear  modulus  relaxation  function.  Again  using 
Table  B-6,  the  shear  modulus  is  given  as  G(t)  = 1/2  E ( t ) / f 1 + v(t)].  As 
discussed  in  step  4,  if  v(t)  is  constant,  G(t)  is  known  directly  from 
E(t) . Otherwise,  G(t)  must  be  determined  by  the  best-fit  procedure 
using  the  data  points:  Gft^)  = 1/2  E(ti>/[1  + v(t^)]. 

The  above  steps  illustrate  one  set  of  transformations  for  a partic- 
ular type  of  e>’  erimental  data.  Other  types  of  transformations  with 
other  types  of  experimental  data  can  be  treated  in  a similar  manner. 
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Appendix  C 


A PERSPECTIVE  OF  VISCOPLASTIC  MODELS 


Throughout;  chis  Appendix,  the  term  "viscoplastic"  is  used  as  a 
general  term  to  embrace  all  constitutive  models  that  contain  elements  of 
elasticity,  viscosity,  and  plasticity.  In  the  main  text  the  discussion 
focused  on  a particular  viscoplastic  model  called  "viscoelast ic-p last ic . " 

The  objective  of  this  Appendix  Is  to  present  a broad  view  of 
viscoplastic  models  and  demonstrate  the  characteristics  of  different 
models.  Concepts  and  imagery  are  emphasized  as  opposed  to  detailed 
equations  and  solution  algorithms. 

To  begin  with,  we  recognize  that  viscoplasticity  is  a phenomo logical 
approach  for  characterizing  the  behavior  of  materials.  That  is,  math- 
ematical models  are  sought  which  can  replicate  the  observed  performance 
of  materials.  Yet,  the  mathematical  models  must  conform  to  certain 
thermodynamic  restrictions. 

Three  mathematical  models  or  constitutive  theories  that  satisfy 
thermodynamic  restrictions  and  are  well  accepted  are:  elasticity, 
viscoelasticity,  and  plasticity.  These  separate  theories  can  (1)  be 
directly  combined  to  provide  viscoplastic  models  (combo-viscoplasticity) 
or  (2)  the  separate  theories  can  be  altered  or  extended  to  provide 
viscoplastic  models  (neo-viscoplast ic ity) . This  is  an  important  dis- 
tinction and  is  the  first  major  division  in  the  hierarchy  ■>•'  v i • • .-plastic 
models  as  shown  in  Figure  C-l  . in  the  right  branch  of  this  figure,  ’■  jf.’ . 
(v)  , and  yL_)  represent  classical  models  for  characterizing  plastic, 
viscous,  and  elastic  response.  Various  combinations  and  arrangements  of 
these  basic  components  define  a particular  combo- v iscoplnst  ic  model, 


Figure  C-l.  Two  classes  of  viscoplastic  models, 


Visual  Model 
(one-dimensional) 


l-AAA— — 


General 

Description 


• Elastic  response 

• Linear  with  oe 

• Time  independent 


• Viscous  response 

• Linear  with  ov 

• Time  dependent 
(no  initial  response) 


• Plastic  response 

• Nonlinear  with  Op, 
dependin  , on  hardening 
function 

• Time  independent 


Figure  C-2.  Components  of  combo-viscoplastic  model. 


such  as  the  viscoelastic-plastic  model  presented  in  the  main  text. 

Other  combo-viscoplastic  models  will  be  illustrated  in  subsequent 
paragraphs . 

The  properties  of  the  classical  models,  (?)  , (?) , and  (?) , are 
illustrated  in  Figure  C-2.  The  identification  of  each  component  may  be 
as  simple  or  complicated  as  needed.  For  example,  (?)  may  be  isotropic 
(two  parameters)  or  anistropic  (three  to  twenty-one  parameters)  . Also, 
the  component  (?)  may  be  isotropic  or  anistropic,  and  each  relaxation 
function  may  have  as  many  parameters  as  required.  Lastly,  (?)  may 
represent  any  classical  plasticity  model,  ranging  from  a simple  non- 
hardening Von  Mises  yield  criterion  (one  parameter)  to  complicated 
"capped"  yield  functions  with  universal  hardening  (e.g.,  twenty-seven 
parameters) . 

The  important  distinction  for  combo-viscoplasticity  is  LhaL  each 
component,  (?)  , (?)  , and  (?) , remains  within  its  classical  limitations. 
That  is,  the  viscous  component,  (?) , deforms  with  time  under  constant 
load  and  responds  in  proportion  to  the  load  magnitude  acting  on  (?) . 

The  plastic  component,  (?)  , deforms  as  a function  of  the  current  load 
magnitude  acting  on  (?)  . There  is  no  explicit  reference  to  time  in  the 
(?)  model  except  implicity  through  the  load. 

Contrast  the  above  to  neo-viscoplasticity  (left  branch  of  Figure 
C-l),  where  the  component  ^?)  may  represent  a modified  , lasticity 
formulation  with  a flow  rule  explicitly  dependent  on  time  (see  for 
example,  Green  and  Naghdi  [25],  Zienkiewicz  and  Cormeau  (5),  and  Perzyna 
14].  Alternatively,  v^VP)  may  represent  a modified  viscoelastic  formula- 
tion, which  accounts  for  plastic-type  responses  by  the  assumption  of 
intrinsic  time.  Effectively,  this  means  real  time  in  portions  of  the 
viscoelastic  model  is  replaced  with  intrinsic  time,  which  is  actually  a 
measure  of  strain.  This  method  is  called  the  "endochronic  theory"  (see 
Bazant  ( 26]  . 

The  relative  merits  of  the  various  neo-viscoplastic  models  will  not 
he  argued  here.  However,  it's  the  author's  contention  that  combo- 
viscoplastic  models  are  preferable  to  the  neo-viscoplastic  models  for 
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che  following  reasons:  (1)  combo  models  provide  a rational  approach  for 
synthesizing  a constitutive  model,  (2)  combo  models  are  formed  from 
well-accepted  constitutive  theories  that  satisfy  thermodynamic  restric- 
tions, (3)  parameter  identification  for  combo  models  is  relatively  easy, 
(A)  pre-existing  elastic,  plastic,  and  viscoelastic  material  information 
can  be  used  directly,  and  (5)  virtually  any  typo  of  observed  material 
behavior  can  be  replicated  by  some  combo-viscoplastic  model. 

There  is  no  limit  to  the  number  of  component  elements  that  can  be 
used  to  synthesize  a combo-viscoplastic  model.  However,  five-element 
models  are  probably  a practical  upper  limit.  As  an  example,  consider 
the  five-element  model  shown  at  the  top  of  Figure  C-3,  along  with  the 
example  four-  and  three-element  models.  Here,  the  models  are  restricted 
to  those  containing  at  least  (?) , (v) , and  (e)  and  have  an  initial 
elastic  response.  (Note:  the  models  shown  are  all  named  using  hyphen 
(-)  to  denote  a series  connection  between  two  elements  and  a slash  (/) 
to  denote  a parallel  coupling.  However,  for  expressing  (v)-(e)  the 
hyphen  is  omitted,  i.e.,  viscoelastic  = visco-elastic,  and  for  (?)-(e) 
we  write  elas toplas t ic  = elastic-plastic.) 

To  appreciate  the  different  behavior  of  these  models,  compare  the 
nature  of  the  viscoelastic-plastic  model  (c)  with  the  elastic- visco/ 
plastic  model  (d)  as  they  respond  (deform)  for  a suddenly  applied  stress. 
If  the  stress  is  less  than  the  elastic  limit  (a  •:  c ),  model  (c)  behaves 
v iscoelastically  and  model  (d)  elastically,  as  shown  in  Figure  C-4.  If 
the  applied  stress  is  greater  than  the  elastic  limit  (::  •»  ) and  the 

hardening  function  t 0,  then  model  (c)  responds  with  an  instantaneous 
clastic  plus  plastic  strain  followed  by  a viscoelastic  creep.  On  the 


o'het  hand,  model  (d)  responds  with  an  instantaneous  elastic  response 
(regardless  of  stress  level)  and  progresses  with  plastic  strain  at  rate 
dependent  on  (v)  . 

As  another  example,  consider  the  tour-element  model  (b) . This 
model  is  capable  of  representing  primary,  secondary,  and  tertiary  creep 
ranges,  as  shown  in  Figure  C-b.  Materials,  such  as  soils,  plastics,  and 
sea-lee,  exhibit  this  type  ol  behavior. 
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Figure  03.  A five-element  combo-viscoplastic  element 
and  its  familv. 
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SYNTHESIZING  COMBO-VISCOPLASTIC  MODELS 


To  synthesize  a viscoplastic  model  (i.e.,  stress-strain  relation- 
ship) from  a particular  configuration  of  (e)  , (v) , and  (?)  requires 
following  a few  simple  rules.  To  wit,  elements  in  series  have  the  same 
stress  and  deform  separately,  whereas  elements  in  parallel  have  the  same 
strain  and  split  the  net  stress  between  them.  In  the  main  text,  these 
rules  were  used  to  synthesize  the  series  model,  viscoelastic-plastic. 
Here  the  synthesis  procedure  is  demonstrated  for  the  parallel  assembly, 
viscoelastic/elastoplastic,  shown  in  Figure  C-6. 

Using  the  notation  in  the  figure,  we  have: 


a +o 
~ep  ~ve 


e + e 
-ep  ~ve 


(C-l) 


(C-2) 


where  a and  e are  the  stress  and  strain  vectors  in  the  elastoplastic 
—ep  -ep 

subassembly,  and  o and  e are  the  stress  and  strain  vectors  in  the 
viscoelastic  subassembly.  Both  series'  subassemblies  were  developed 
fully  in  the  main  text.  Thus,  from  Equation  2-33,  the  elastoplastic 
stress-strain  relationship  is: 


D )Ac 

" P - 


ep 


(C— 3 ) 


where  Dj  is  the  elastic  matrix  representing  , and  is  the  plas- 
ticity matrix,  which  is  a function  of  o (i.e.,  not  total  stress  o) . 

-ep 

Similarly,  from  Equation  3-22,  the  viscoelastic  stress-strain  relation- 
ship is: 
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(D,  + D )Ac  + a 

»1  «.v  ~ve  ~v 


where  is  the  elastic  matrix  representing  (E^) , Dy  is  the  viscoelastic 
matrix  dependent  on  time  step,  and  £y  is  the  stress  history  matrix. 

Adding  Equations  C-3  and  C-4  and  using  the  relationship,  in  Equa- 
tions C-l  and  C-2 , the  final  stress-strain  relationship  is: 


Ao  = (D  + D - <D  >)  Ac  - 5 (C-5: 

0+  V m P 

In  the  above,  Dg  = ^ + p2  is  the  total  elastic  matrix.  If  Dj.  = 0,  the 
model  degenerates  to  an  elastoplastic  model,  or  if  D,  = 0,  the  model  is 
viscoelastic.  Any  other  linear  combination  of  pi  and  p2  provides  a 
viscoplastic  model. 

Other  combo-viscoelastic  models  can  be  developed  by  following  the 
procedures  outlined  above. 


Figure  C-6.  Viscoelastic/elastoplastic  model. 
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NAS  CO.  Guani.m.inio  Bav  Cuba.  Code  114.  Alameda  CA:  Code  l!0  I Fac.  Plan  HR  MOKi.  ( ode  18’.  Jackson » cilc  IT  . 
Code  HOIK).  Brunswick  M I..  Code  A244  iG.  'I  iaski.  Pom:  Mugu  CA ; Code  T).  Atlanta.  Marietta  ( I A:  Du  lnl  Div  . 
Bermuda;  I NS  Uuthhol/.  Pensacola.  11.1  v.nl  Chief.  Petty  Ofir  PW  Self  Help  Div.  Bccvillc  I \.  I”A  -I 
MagtiilC).  C orpio  t hnsli  IX.  TWI)  Maim.  Div  . New  Orleans.  Belle  Chevse  LA:  PW  I),  W illow  (ire1.-  I ’ X PWt ) 
iM  F.llmli)  Los  Alamitos  CA.  PV’.t)  Belle  V havse.  I A.  I’W  ( ) t hase  I icld  Becvillc.  TX . PW O Key  ’ > 'i  f I . PW 1 1 
W hiting  Eld.  Milton  IT.:  PWT.-  ! 5all.iv  1 X ; PW  ().  Glenview  1 1 . PW  ( ).  Kmgw die  I X.  PW  ( i.  Mn.iui.ti  S.m  1 tie-go 
CA;  SCE  L-int  Heel  Norfolk  VA.  SCI  Norfolk.  V A.  St  I Barbers  Point  HI 
N A I L R ('.SEARCH  COL  NCI  I . N.iv.d  Studies  Board.  W ashing  ton  DC 
NA  TPARACTH  I I I LS  I RAN  PW  F.ngi.  FT  Centro  CA 
NAVACI  PWO,  London  UK 
NAVAF.ROSPRI  GMI  IK  I N St  F.  Penv.,eo|a  I I 

NAVAL  FACILITY  PWO  Barbados.  PWO.  Brawdy  Wales  UK.  PWO.  C ape  Hailcr.is.  Buvlon  NC 
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NAVCOASTSYSLAB  CO.  Panama  City  FL;  Code  423  (D.  Good),  Panama  City  FL;  Code  715  (J.  Quirk)  Panama  City. 
FL:  Library  Panama  City,  FL 

NAVCOMMAREAMSTRSTA  PWO.  Norfolk  VA;  PWO.  Wahiawa  HI;  SCE  Unit  I Naples  Italy 
NAVCOMMSTA  Code  401  Nea  Makri.  Greece;  PWO,  Adak  AK;  PWO,  Exmoulh,  Australia 
NAVCONSTRACEN  COtCDR  C.L.  Neugent),  Port  Hueneme,  CA 
NAVEDTRAPRODEVCEN  Tech.  Library 
NAVEDUTRACEN  Engr  Dept  (Code  42)  Newport,  RI 
NAVEODFAC  Code  605,  Indian  Head  MD 

NAVFACENGCOM  Code  043  Alexandria.  VA;  Code  044  Alexandria.  VA;  Code  0451  Alexandria.  VA;  Code  0453  (D. 
Potter)  Alexandria.  VA;  Code  0454B  Alexandria.  Va;  Code  04B5  Alexandria,  VA;  Code  101  Alexandria.  VA,  Code 
10133  (J.  Leimanis)  Alexandria.  VA;  Code  1023  (T.  Stevens)  Alexandria.  VA;  Code  104  Alexandria.  VA;  Code  2014 
(Mr.  Taam).  Pearl  Harbor  HI;  Morrison  Yap,  Caroline  Is.;  P W Brewer  Alexandria.  VA;  PC-22  (E.  Spencer) 
Alexandria.  VA;  PL-2  Ponce  P R.  Alexandria.  VA 

NAVFACENGCOM  - CHF.S  DIV.  Code  101  Wash.  DC,  Code  402  (R.  Morony)  Wash.  DC;  Code  405  Wash,  DC;  Code 
FPO-I  (C  Bodey)  Wash.  DC;  Code  FPO-I  (Ottsen)  Wash,  DC,  Code  FPO-ISP  (Dr.  Lewis)  Wash,  DC;  Code 
FPO  1SP13  (T  F Sullivan)  Wash.  DC;  Code  FPO-IPI2(Mr.  Scola).  Washington  DC;  Scheessele,  Code  402.  Wash. 
DC 

NAVFACENGCOM  - LANT  DIV  ; Eur.  BR  Deputy  Dir.  Naples  Italy;  RDT&ELO09P2.  Norfolk  VA 
NAVFACENGCOM  - NORTH  DIV  (Boretsky)  Philadelphia.  PA;  Code09P(LCDR  A J.  Stewart);  Code  1028. 

RDT&F.LO.  Philadelphia  PA;  Design  Div.  (R.  Masino).  Philadelphia  PA;  ROICC.  Contracts,  Crane  IN 
NAVFACENGCOM  - PAC  DIV  Code  402.  RDT&E.  Pearl  Harbor  HI;  Commander.  Pearl  Harbor.  HI 
NAVFACENGCOM  - SOUTH  DIV  Code  90.  RDT&ELO.  Charleston  SC;  Dir.,  New  Orleans  LA 
NAVFACENGCOM  - WEST  DIV.  Code  04B;  O9P/20;  RDT&ELO  Code  2011  San  Bruno.  CA 
NAVFACENGCOM  CONTRACT  AROICC.  Point  Mugu  CA;  AROICC.  Quantico,  VA;  Eng  Div  dir.  Southwest  Pac. 
Manila.  PI:  OICC.  Southwest  Pac.  Manila.  PI;  OICC/ROICC.  Balboa  Canal  Zone:  ROICC  (Ervin;  Puget  Sound 
Nasal  Shipyard.  Bremerton.  W'A.  ROICC  AF  Guam;  ROICC  LANT  DIV  . Norfolk  VA:  ROICC  Ofl  Point  Mugu. 
CA.  ROICC.  Diego  Garcia  Island:  ROICC.  Keflavtk.  Iceland;  ROICC,  Pacific.  San  Bruno  CA 
NAVHOSP  LT  R.  Elsbernd.  Puerto  Rico 
NAVMAG  SCF.  Guam 
NAVMIRO  OIC.  Philadelphia  PA 

NAVOCEANSYSCEN  Code  409  (D.  G Moore).  San  Diego  CA;  Code  531  l(T)  (E  Hamilton)  San  Diego  CA;  Code 
6565  (Tech.  Lib.),  San  Diego  CA:  Research  Lib,.  San  Diego  CA 
NAVORDSTA  PWO,  Louisville  KY 
NAVPETOFF  Code  30.  Alexandria  VA 
NAVPGSCOL  Code  61 WL  (O  Wilson)  Monterey  CA 

NAVPHIBASE  CO.  ACB  2 Norfolk.  VA  Code  S3T.  Norfolk  VA;  Harbor  Clearance  Unit  Two.  Little  Creek,  VA; 

OIC.  UCT  ONE  Norfolk.  Va 

NAVREGMEDCEN  Code  3041,  Memphis.  Mtlltnglon  TN;  SCE  (D  Kavt).  SCEtLCDR  B.  E.  Thurston).  San  Diego 
CA;  SCF.  Camp  Pendleton  CA 

N AVSCOl.CECOFF  C35  Port  Hueneme.  CA;  C44A  (R  Chittenden).  Port  Hueneme  CA;  CO,  Code  C4aA  Port 
Hueneme.  CA 

NAVSEASYSCOM  Code  OOC  (LT  R.  MacDougal).  Washington  DC 
NAVSEC  Code  6034  (Library).  Washington  DC 
NAVSECGRUACT  PWO.  Torn  Sta.  Okinawa 
NAVSHIPREPFAC  Library.  Guam.  SCE  Subic  Bay 

NAVSH1PYI3.  Code  202  4,  Long  Beach  CA.  Code  202.5  ( Library ) Fhiget  Sound.  Bremerton  W'A:  Code  4UO.  Puget 
Sound.  Code  404  (LT  J Rtccto).  Norfolk.  Portsmouth  VA;  Code  410.  Mare  Is..  VallejoCA.  Code  440  Portsmouth 
NH;  Code  440.  Norfolk;  Code  440.  Puget  Sound.  Bremerton  W'A;  Code  440.4,  Charleston  SC;  L.D.  Vivian. 

Library.  Portsmouth  NH;  PWO.  Mare  Is  ; Tech  Library.  Vallejo.  CA 
NAVSTA  CO  Naval  Station.  Mayport  FL.  CO  Roosevelt  Roads  P R.  Puerto  Rico;  Engr.  Dir.,  Rota  Spain:  Maint.  Di\ 
Dir- Code  531 . Rodman  Canal  Zone;  PW  D (l.TJG.P  M Moiolenichi.  Puerto  Rico;  PWO  Midway  Island;  PWO. 
Keflavik  Iceland;  PWO.  Mayport  FL;  ROICC.  Rota  Spain;  SCE.  Guam;  SCE,  Subic  Bay.  R P : Utilities  Engr  Off. 
tLTJG  A.S.  Ruchiei,  Rota  Spam 
NAVSU BASF.  LTJG  D.W  Peck.  Groton.  CT 

NAV'SUPP.ACT  CO.  Seattle  W'A.  Code  413,  Seattle  W'A;  Engl.  Div  <F  Mollica).  Naples  Italy;  LTJG  McGarrah. 
Vallejo  CA 

NAVSU RFWPNCFN  PWO.  White  Oak,  Silver  Spring.  MD 
NAVTEGHTRACEN  SCE.  Pensacola  FL 

NAVW  PNUF.N  Code  2636  |W.  Bonner).  China  Lake  CA;  PWO  (Code  26).  China  Lake  CA;  ROICC  (Code  702).  Chin  i 
Lake  CA 

NAVW  PNSTA  EARLL  FNS  G.A.  Lowry . Fallbrook  CA;  Maint.  Control  Dir..  Yorktown  VA:  PW'  Office  (Code  09CI I 
Yorktown.  VA;  Security  Offr,  Colts  Neck  NJ 
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NAVWPNSUPPCEN  Code  09  (Boenmghausen)  Crane  IN 
NCBU  405  OIC.  San  Diego,  CA 

NCBC  CEL  (CAPT  N.  W.  Petersen),  Port  Hueneme.  CA;  CEL  AOIC  Port  Hueneme  CA;  Code  10  Davisville.  Rl; 
Code  155,  Port  Hueneme  CA;  Code  156,  Port  Hueneme,  CA;  PW  Engrg.  Gulfport  MS;  PWO  (Code  80)  Port 
Hueneme,  CA 

NCBU  411  OIC,  Norfolk  VA 
NCR  20.  Commander 

NMCB  I33(ENS  T.W.  Nielsen);  5,  Operations  Dept.;  74.  CO;  Forty,  CO;  THREE,  Operations  Off. 

NORDA  Code  440  (Ocean  Rsch  Off)  Bay  St.  Louis  MS 

NRL  Code  8400  (J.  Walsh).  Washington  DC;  Code  8441  (R.A.  Skop).  Washington  DC;  Rosenthal,  Code  8440.  Wash. 
DC 

NSC  Code  54.1  (Wynne).  Norfolk  VA 
NSD  SCE.  Subic  Bay.R.P. 

NTC  Code  54  (F.NS  P.  G.  Jackel).  Orlando  FL;  Commander  Orlando.  FL 

NUSC  Code  131  New  London.  CT;  Code  EAI23  (R.S.  Munn),  New  London  CT , Code  TAI3I  (G.  De  la  Cruz),  New 
London  CT 

ONR  Code  700F  Arlington  VA;  Dr.  A.  Laufer,  Pasadena  CA 
PHIBCB  I P&E,  Coronado.  CA 

PMTC’  Code  4253-3,  Point  Mugu,  CA,  Pat.  Counsel,  Point  Mugu  CA 

PWC  ENS  J.F..  Surash.  Pearl  Harbor  HI,  ACE  Office  (LTJG  St.  Germain)  Norfolk  VA;  CO  Norfolk.  VA.  CO.  Great 
Lakes  IL;  Code  1 16  (LTJG  A.  Eckharl)  Great  Lakes.  1L;  Code  120.  Oakland  CA;  Code  120C  (Library)  San  Diego. 
CA.  Code  1 28.  Guam.  Code  200.  Great  Lakes  IL;  Code  200,  Oakland  C A;  Code  220  Oakland.  CA;  Code  220  I , 
Norfolk  VA;  Code  30C  (Boettcher)  San  Diego.  CA;  Code  680,  San  Diego  CA.  Library.  Subic  Bay,  R P ; OIC 
CBL-405.  San  Dtego  CA.  XO  Oakland,  CA 
SPCC  Code  122B,  Mechamcsburg.  PA;  PW'O  (Code  120)  Mechanicsburg  PA 
U S MERCHANT  M.4RINL  ACADEMY  Kings  Point.  NY  (Reprint  Custodian) 

US  GEOLOGICAL  SURVEY  Off.  Marine  Geology.  Piteleki,  Reston  VA 

USAF  SCHOOL  OF  AEROSPACE  MEDICINE  Hyperbaric  Medicine  Div.  Brooks  AFB.  TX 

USCO  (G-ECV)  W ashington  Dc;  (G-ECV/61 ) (Burkhart)  Washington.  DC.  G-EOE-4/61  (T.  Dowd).  Washington  DC 

USCG  ACADEMY  LT  N Stramandi.  New  London  CT 

USCG  R&D  CENTER  D.  Moiherway,  Groton  CT;  Tech.  Dtr.  Groton,  CT 

USNA  Ch.  Mech.  Engr.  Dept  Annapolis  MD;  Ocean  Sys.  Eng  Dept  (Dr.  Monney)  Annapolis.  MD;  PWD  Engr.  Div. 

(C.  Bradford)  Annapolis  MD;  PWO  Annapolis  MD 
AMERICAN  CONCRETE  INSTITUTE  Detroit  Ml  (Library  ) 

CALIF  DEPT  OF  NAVIGATION  & OCEAN  DEV.  Sacramento.  CA  (G.  Armstrong) 

CALIFORNIA  STATE  UNIVERSITY  LONG  BEACH.  CA  (CHELAPAl !);  LONG  BEACH.  CA  (YEN) 
COLORADO  STATE  UNIV..  FOOTHILL  CAMPUS  Fort  Collins  (Nelson) 

CORNELL  UNIVERSITY  Ithaca  NY  (Serials  Depi.  Engr  Lib.) 

DAMES  & MOORE  LIBRARY  LOS  ANGELES,  CA 

DUKE  UNIV  MEDICAL  CENTER  B.  Muga.  Durham  NC;  DURHAM.  NC  (VESIC) 

FLORIDA  ATLANTIC  UNIVERSITY  BOCA  RATON.  FI.  (MC  ALLISTER).  Boca  Raton  FL  (Ocean  Engr  Dept  . C. 
l.in) 

FLORIDA  ATLANTIC  UNIVERSITY  Boca  Ralon  FL(W.  Tessin) 

FLORIDA  TECHNOLOGICAL  UNIVERSITY  ORLANDO.  FL  (HARTMAN) 

GEORGIA  INSTITUTE  OF  TECHNOLOGY  Atlanta  GA  (School  of  Civil  Engr,.  Kahn):  Atlanta  GA  (B  Mazanti) 
IOW  A STATE  UNIVERSITY  Ames  IA  (CE  Depi,  Handy) 

VIRGINIA  INST  OF  MARINE  SCI.  Gloucester  Point  VA  (Library  ) 

LEHIGH  UNIVERSITY  BETHLEHEM,  PA  (MARINE  GEOTECHNICAL  LAB..  RICHARDS):  Bethlehem  PA 
( Fritz  Engr.  Lab  No.  13.  Bcedle).  Bethlehem  PA  (l.tnderman  l.ib.  No. 30.  Flecksletner) 

LIBRARY  OF  CONGRESS  WASHINGTON.  DC  (SCIENCES  & TECH  DIV) 

MICHIGAN  TECHNOLOGICAL  UNIVERSITY  Houghlon.  Ml  (Haas) 

Mi  l Cambridge  MA:  Cambridge  MA  (Rm  10-500.  Tech.  Reports.  Engr.  Lib  );  Cambridge  MA  (Whitman) 

NEW  MEXICO  SOLAR  ENERGY  INST  Dr  Zwibel  Las  Cruces  NM 

NORTHWESTERN  UNIV  Z P Bazant  Evanston  IL 

NY  CITY  COMMUNITY  COLLEGE  BROOKLYN.  NY  (LIBRARY) 

UNIV.  NOTRE  DAME  Kutona.  Notre  Dame,  IN 

OREGON  STATE  UNIVERSITY  (CE  Dept  Grace)  Corvallis.  OR;  CORVALLIS.  OR  (CE  DFPT.  BELLI; 

CORVALLIS.  OR  (CE  DEPT.  HICKSi;  Corvalis  OR  (School  of  Oceanography) 

PENNS'!  LVANIA  STATE  UNIVERSITY  UNIVERSITY  PARK.  PA  (GOTOLSKI) 

PU  RDUE  UNIVERSITY  Lafayette  IN  (Leonards),  Lafayette.  IN  ( Altschaeffll;  Lafayette.  IN  (CE  Engr  l ib) 

SAN  DIEGO  STATE  UNIV.  I.  Nooruny  San  Diego.  CA;  Dr.  Krishnamoorthy . San  Diego  CA 


SEATTLE.  U Prof  Sel.wjegler  Seattle  WA 

SOUTHWEST  RSCH  INST  King.  San  Aniomo,  TX;  R.  DeHart.  San  Antonio  TX 
STANFORD  UNIVERSITY  Engr  Lib.  Stanford  C A;  Stanford  CA  (Genet 
STATE  UN1V.  OF  NEW  YORK  Buffalo.  NY 

TEXAS  A&M  UNIVERSITY  COLLEGE  STATION.  TX  (CF.  DEPT);  College  Station  TX  tCF.  Dept.  Herbich) 
UNIVERSITY  OF  CALIFORNIA  BERKELEY.  CA  (CF.  DEPT.  GKRWICK);  BERKELEY,  CA  (CE  DEPT. 
MITCHELL);  Berkeley  CA  (B.  Bresler);  Berkeley  CA  (Dept  of  Naval  Arch  ).  Berkeley  CA  (R.  Williamson); 
DAVIS.  CA  (CE  DF.PT.  TAYLOR);  LIVERMORE.  CA  (LAWRENCE  LIVERMORE  LAB.  TOKARZ);  La  Jolla 
CA  t.Acq.  Dept.  Lib  C-075.A).  Los  Angeles  C A (F.ngr  I,  K.  Lee).  M Duncan.  Berkelev  CA;  SAN  DIEGO,  CA,  l.A 
JOLLA.  CA  (SF.ROCKI) 

UNIVERSITY  OF  DELAW  ARE  Newark.  DE  iDept  of  Civil  Engineering,  Chcsson) 

UNIVERSITY  OF  HAW  AII  Dr  Chiu  Honolulu.  HI  HONOLULU.  HI  (SCIENCE  AND  TECH.  DIV.);  Honolulu  HI 
I Ur.  Siilardl 

UNIVERSITY  OF  ILLINOIS  Metz  Ref  Rm.  Lrbana  IL;  URBAN  A.  II.  (DAVISSON);  URBAN  A,  It  (LIBRARY  I: 

URBANA.  IL  (NEW  ARK);  Urban;.  ILiCF  Dept.  W Gamble) 

UNIVERSITY' OF  MASSACHUSETTS  (Heronemus).  Amherst  MACE  Dept 
UNIVERSITY  OF  MICHIGAN  Ann  Arbor  Ml  (Richard 

UNIVERSITY  OF  NEBRASKA -LINCOLN  Lincoln.  NF.  (Ross  Ice  Shelf  Proj.) 

UNIVERSITY  OF  NEW  MEXICO  J Nielson-Engr  Malls  & Civ li  Sys  Div.  Albuquerque  NM 
UNIVERSITY  OF  TEXAS  Inst.  Marine  Sei  i Library  ).  Port  Arkansas  TX 

UNIVERSITY  OF  TEXAS  AT  AUSTIN  AUSTIN.  TX  (THOMPSON i;  Austin  TX  <R  Olson);  Auslin.  TX  (Bieenl 
UNIVERSITY  OF  WASHINGTON  Dept  of  Civil  Fngr  (Di.  Mattock).  Seattle  WA;  SFA1TUK.  W A (MERCHANT): 
SEATTLE.  WA  UK  F AN  ENG  RSCH  LAB.  GRAY);  Seattle  W A (E.  Linger);  Seattle.  W A Transportation. 
Construction  & Geom  Div 

CHS  RESEARCH  CO.  LIBRARY  SAN  MATEO.  CA 

ALFRED  A.  YEF  <1  ASSOC  Honolulu  HI 

AMF.TEK  Offshore  Res  & Fngr  Div 

APPLIED TECH  COUNCIL  R.  Scholl.  Palo  Alio  CA 

AH V 111  GRANT  OLYMPIA.  WA 

ATLANTIC  RICHFIELD  CO  DALLAS.  TX  (SMITH! 

AUSTRALIA  Dept  PW  (A.  Hicks).  Melbourne 
BECHTEL CORP.  SAN  ERANCISCO.CAtPHEl.PS) 

RIT.GIUM  HAEC'ON,  NY.  Gent 
BROWN  A ROOT  Houston  T X (D.  Ward) 

CANADA  Mem  Unis  Newfoundland  iChani.  St  Johns;  Surveyor,  Nenninger  A Chcnevcrt  Inc  . Montreal;  Warnock 
Hersey  Prof  Srv  l td.  La  Sale.  Quebec 
CE  BRAUN  CO  llu  Bonchet.  Murray  Hill.  NJ 

CHEVRON  OIL  FIELD  RESEARCH  CO  LA  HABRA,  CA  (BROOKS) 

CONCRETE  TECHNOLOGY  CORP  TACOMA.  WA  t ANDERSON) 

CONRAD  ASSOC.  Van  Nuys  C.A  l.A.  I.uisnni) 

DR  A VO  CORP  Pittsburgh  PA  (Giannino!;  Pillsburgh  PA  (Wrightl 

NORW  AY  OFT  NORSK!:  VERIT  AS  tl.ibiary  ).  Oslo 

KVAI.U AVION  ASSOC  INC  KING  OF  PRUSSIA.  PA  (IT  DELE) 

FORD.  BACON  & DAVIS.  INC  New  York  ( Library  ) 

FRANCE  Dr  Outer  Ire.  Boulogne.  1 . Pliskin.  Pans.  P.  Jensen.  Boulogne 
CLOT T CHNICAL  ENGINEERS  INC  Winchester.  MA  iPaulding) 

(iLIDDFNCO.  ST  RONGSVII  LI  . OH  (RSCH  LIB) 

GRUMMAN  AEROSPACE  CORP.  Bethpagc  NY  (Tech.  Into.  Ctrl 
HM  I Y A ALDRICH.  INC.  Cambridge  MA  t Aldrieh.  Jr  ) 

H( )N  E.Y  YY  El  I , INC  Minneapolis  M N (Residential  F.ngr  Lib  ) 

HUGHES  AIRCRAFT  Culvci  City  CA  (Tech.  Doe  Ctrl 

IT  Yl  Y M C.nroni.  Milan;  Sergio  Tuiion:  Mil.no;  Torino  (F.  Levi) 

MAKAI  OCEAN  LNGRNG  INC  K.ulu.i.  HI 
JAMES  CO.  K Girdles.  Oilando  El 

1 AMON  I IXJHEKI  Y Gl  OEOOICAl  OBSI.RY  Palisades  NY  i McCoy i;  1’alisai.cs  NY  (Selwyni 
LOCK  H KIT!  Ml  SSI I I S A SPACE  CO.  INC.  Mgt  Nav  :tl  Arch  A Mar  Eng  Sunny  vale.  CA.  Sunnyvale  CA 
iRynewic/l.  Sunnyvale.  CA  (Phillips) 

1 OCKHI  I DOC  I AN  l ABOKATORY  San  Diego  CA(F  Simpson) 

MARAT  HON  Oil  CO  Houston  IX  (C  Seay  I 

MC'  C l .1.1  I AND  ENGINEERS  INC  Houston  IX  t».  MeClellmd) 

MC  DONNE  I AIRCRAFT  CO.  Depl  .Mil  ,K  H.  Faymani.  St  l unis  Mo 
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MEDALL  & ASSOC.  INC  J.T  GAEFEY  ll  SANTA  ANA.  CA 
MEXICO  R.  Cardenas 

MOBIL  PIPE  LINE  CO.  DALLAS.  TX  MGR  OF  ENGR  <NOACK> 

MUESER.  RUTLEDGE.  WENTWORTH  AND  JOHNSTON  NEW  YORK 'RICHARDS) 

NEW  ZEALAND  New  Zealand  Concrete  Research  Assoc.  (Librarian).  Porirua 
NEWPORT  NEWS  SHIPBLDG  & DRY  DOCK  CO  Newport  News  V A (Tech.  I ih. ) 

NORW  AY  DET  NORSKE  VERITAS  tRoren)  Oslo:  I.  Foss.  Oslo;  J.  Creed,  Ski:  Norwegian  Tech  Univ  f Brandtzaeg), 
Trondheim 

OFFSHORE  DEVELOPMENT  ENG  INC.  BERKELEY.  CA 
PACIFIC  MARINE  TECHNOLOGY  LONG  BEACH.  CA  iW  AGNFR) 

PORTLAND  CEMENT  ASSOC.  SKOKIE.  ILtCORELY):  Skokie  IL  (Rvh&  Dev  L.ah.  Lib.) 

PRESCON  CORP  TOWSON.  MD  (KF.LLER) 

RAND  CORP  Santa  Monica  CA  (A  Laupa) 

RAYMOND  INTERNATIONAL  INC.  E Colie  Soil  Tech  Depl.  Pennsauken.  NJ 
RIVERSIDE  CEMENT  CO  Riverside  CA  (W.  Smilh) 

SAND1A  LABORATORIES  Library  l)i\  ..  Livermore  CA 
SCHUPACK  ASSOC  SO  NORW  ALK.  CT  tSCHUPAC  Ki 
SEATECH  CORP.  MIAMI.  FL  (PERONI) 

SHELL  DEVELOPMENT  CO.  Houston  TX(E.  Doyle) 

SHELL  OIL  CO.  HOUSTON.  TX  (MARSHALL) 

SOUTH  AMERICA  N.  Nouel.  Valencia.  Venezuela 
SWEDEN  GcoTech  Inst;  VBB  (l  ibrary).  Stockholm 
TIDEWATER  CONSTR  CO  Norfolk  VA  iFowlerl 

1TCW  SYSTEMS  CLEVELAND.  OH  tENC,  LIB  ):  REDONIM)  BEACH.  CA  (DAI) 

UNITED  KINGDOM  Cement  & Concrete  Assoc  V^esham  Springs.  Slough  Bucks:  Cement  & Concrete  Assoc.  (Lit. 
Em.  Bucks;  D.  Lee.  London;  D New.  G.  Maunsell  & Partners,  London;  Library . Br  ol;  Shaw  A Hatton  iE 
Hansen).  London:  Taylor.  Woodrow  Constr  (0I4P).  Southall.  Middlesex;  Univ  ol  Bristol  tR.  Morgan i.  Ruslol 
WATT  BRIAN  ASSOC  INC.  Houston.  TX 

W F.STINGHOUSF  ELECTRIC  CORP.  Anr.apolis  MD  (Oceanic  Div  Lib.  Bryan):  Library.  Pittsburgh  PA 
WISS.JANNKY.  EL.STNER.  & ASSOC  Northbrook.  1L  (J.  Hanson) 

WOODWARIK'LYDF.  CONSULTANTS  (A  Harrigan)  San  Francisco;  PLYMOUTH  MEETING  PA  (CROSS.  Ill) 
AL  SMOOTS  l.os  Angeles.  CA 
BARA.  JOHN  P Lakewood.  CO 
BROWN.  ROBERT  University.  Al. 

BULLOCK  La  Canada 
E.  HEL'ZK  Boulder  CO 
CAP!  MURPHY  Sunnyvale,  CA 
GREG  PAC.F.  EUGENE.  OR 
R E BESIF.R  Old  Saybrook  CT 
T V. . MF.RMEL  Washington  DC 
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